Monitoring rangelands by identifying the departure of contemporary conditions from long-term ecological potential allows for the disentanglement of natural biophysical gradients driving change from changes associated with land uses and other disturbance types. We developed maps of ecological potential (EP) for shrub, sagebrush (Artemisia spp.), perennial herbaceous, litter, and bare ground fractional cover in Wyoming, USA. EP maps correspond to the potential natural vegetation cover expected by environmental conditions in the absence of anthropogenic and natural disturbance as represented by the greenest and least disturbed period of the Landsat archive. EP was predicted using regression tree models with inputs of soil maps and spectral data associated with the 75th percentile of the Normalized Difference Vegetation Index in the Landsat archive. We trained our EP models with 2015 component cover maps on ecologically intact sites with relatively lower bare ground than expected. We generated departure of vegetation cover by comparing the EP and 2015 fractional cover. The departures represent land cover change from potential land cover and/or within-state changes in 2015. Next, we converted EP and 2015 fractional cover maps into thematic land cover and evaluated departure to determine if it was great enough to result in land cover change. The 2015 conditions showed reduced shrub, sagebrush, litter, and perennial herbaceous cover and increased bare ground relative to EP. Known disturbances, such as energy development, fires, and vegetation treatments, are clearly visible on the departure maps, but not on EP component maps. The most frequent departure from EP land cover was shrubland conversion to grassland. Land cover departures can be explained only in small part by known disturbance, and instead are ostensibly related to climate and land management practices. These drivers result in land cover departures that broadened the ecotone between shrubland and grassland relative to EP.
Introduction
Ecosystems naturally vary across space and time due to biophysical properties such as soils, topography, hydrology, and climate. Ecosystems also vary spatially and temporarily due to anthropogenic drivers, such as introduction of invasive species, increased disturbances to the land, and alterations of ecosystem processes with short- and long-term effects on the ecological potential. Ecological potential (sometimes referred to as site potential) is the possible natural vegetation cover and type expected by environmental conditions in the absence of anthropogenic and natural disturbance (Hickler et al., 2011; Gu and Wylie 2010; Rigge et al., 2019). Ecological potential (EP) data are critical then, to serve as a baseline for comparison to current conditions. Some of the many factors that alter ecosystem productivity include soil erosion, compaction, nutrient depletion, loss of soil organic carbon, salinization, climate, management, invasive species, fires, and flooding (Van der Esch et al., 2017). Globally, Van der Esch et al. (2017) estimated that 23% of global land area has productivity lower than its potential natural state. There is wide disagreement, however, in estimates of total degraded land area, which vary from 1 to over 6 billion ha globally (Gibbs and Salmon, 2015). Moderate or greater loss of biotic integrity (e.g. productivity) from potential conditions was found in 21.3% of evaluated rangelands over the western United States (U.S.) (Herrick et al., 2010). Ecological potential data are useful to determine if an ecosystem can recover to a desired state following disturbance without management intervention (Standish et al., 2014), set restoration goals, and adapt management practices to changing environmental conditions (White and Walker, 1997).
Two methods for estimating the EP include using contemporary data (spatial potential) from undisturbed sites that are analogous to the site in question (White and Walker, 1997) and using historical data (temporal potential) at a given site. Spatial potential can be based on expert opinion, potential natural vegetation, or remote sensing (Herrick et al., 2019). However, both approaches have drawbacks that add cost and uncertainty to land management projects (Herrick et al., 2019). Spatial potential methods often suffer from ambiguity in selecting appropriate reference sites and may not represent dominant conditions. For example, dry sites in tallgrass prairie are the least likely to have been converted to agricultural use, which would result in a scarcity of references in more mesic tallgrass sites and overrepresentation of dry sites. In other cases, undisturbed conditions may not be the EP or the vegetation type may not exist currently (Friedel, 1991). Temporal potential methods often have incomplete data through time and may exclude unmeasured vegetation transitions (e.g., fire) that confound change interpretation (White and Walker, 1997). Temporal potential methods would ideally span multiple generations of dominant plant species. Using a short temporal reference can either ignore or accentuate rare events (e.g., fire or very wet years), which may be a prerequisite to recovery (Friedel, 1991) or diminish the importance of interannual weather variation. The often limited spatial and temporal extent of relevant field data further limits broad assessments of ecosystem condition and resilience.
Remote sensing data are spatially extensive, and depending on the sensor, have a rich and consistent archive of imagery extending back decades, both critical to accurate mapping of ecosystem potential and current conditions. Assessing the contemporary conditions of a rangeland site often occurs through understanding the factors affecting soil stability, hydrologic function, and biotic integrity (i.e., Interpreting Indicators of Rangeland Health: Pyke et al., 2002), an approach that has been widely adopted (Herrick et al., 2019). Key indicators of these functions are bare ground cover, soil resistance to erosion, plant structural groups, annual productivity, litter cover, and presence of invasive plants (Pyke et al., 2002), which can all be inferred directly or indirectly using remote sensing techniques (McCord et al., 2017, Xian et al., 2015). Additional methods for identifying degradation include the detection of opposing Net Primary Productivity (NPP) and precipitation trends (Yengoh et al., 2014, Cowie et al., 2018).
Remote sensing approaches often rely on the Normalized Difference Vegetation Index (NDVI) (Herrick et al., 2019). However, in some cases, degradation processes may increase NDVI (Bai et al., 2008), such as shrub encroachment or expansion of invasive herbaceous grasses (Herrick et al., 2019). Some ecosystems with frequent fire such as the tallgrass prairie could often bioclimactically support forest, as evidenced by the widespread invasion of Juniperous spp. (e.g., Briggs et al., 2005; Wang et al., 2018), In such ecosystems, EP would be subject to the length of the reference period, and its interpretation would require local knowledge. Reeves and Bagget (2014) identify a process to assess the trend and state of vegetation production relative to a potential condition by comparing Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI at neighboring pixels of similar potential communities, structure, and productivity to the site of interest. Rigge et al. (2019) developed an approach that leverages the Landsat archive to produce the 90th percentile NDVI over the 1985–2016 period combined with biophysical data to model long-term EP productivity using regression trees (RT). These data were compared to 2015 component cover maps (Xian et al., 2015; Rigge et al., 2020a) for rangelands across the Great Basin of the western U.S. to produce component cover departure. Areas of departure were associated with fires, land management practices, and other disturbances in a manner readily understood by land managers. Departure data are invaluable as rangeland component cover data are related to ecosystem and landscape structure and function (Xian et al., 2015).
Our primary objective was to produce EP maps of rangeland component cover (shrub, sagebrush [Artemisia spp.], herbaceous, litter, and bare ground) for the state of Wyoming, USA based on the least disturbed and best-growing conditions in the Landsat archive, trained on ecologically intact sites. Best-growing conditions are those commensurate with the 75th percentile of NDVI and concurrent spectral bands as observed in each growing season month in the Landsat archive (1984–2018). Next, we compared these EP component cover maps to actual 2015 component cover maps (Rigge et al., 2020a) to generate differences from potential in units of percent cover. These differences (EP departure) represent departures from potential land cover and/or within-state changes. Our second objective was to evaluate ecosystem and regional patterns in component cover departure, the role of climate, and of 2015 weather conditions to determine the relevance of departure maps in management. Third, we converted EP and 2015 component cover maps into thematic land cover to evaluate land cover change, with the goal of exploring whether EP component cover from the Landsat archive and 2015 component cover can detect shifts in the shrubland/grassland boundary and account for natural and anthropogenic changes on the landscape that occur at multiple spatial and temporal scales. We hypothesized that 1) 2015 component cover will have lower vegetation cover than EP, and therefore 2) the shrubland/grassland boundary will shift.
Study site
Our study area includes Wyoming rangelands, and portions of Montana, South Dakota, Nebraska, Colorado, and Utah, mapped by Rigge et al. (2020a), totaling 201,464 km2 (Fig. 1). Study area average precipitation ranges from 112 to 726 mm, with a mean of 354 (Thornton et al., 2016), and elevation ranges from 956 to 4013 m. Fires are relatively rare in the region, with only 3.8% of the study area burned at least once from 1984–2017, most fires occurred in the north and east. Portions of five U.S. Environmental Protection Agency (EPA) Level III ecoregions are present in the study area (Fig. 1). The Wyoming Basin ecoregion is dominated by stands of sagebrush steppe with a modest herbaceous understory and a semi-arid environment. Isolated mountain ranges and riparian areas interrupt the sagebrush plain, containing a greater abundance of herbaceous cover. The Northwestern Great Plains ecoregion has less shrub overall than the Wyoming Basin, and an overall pattern of declining shrub cover from west-to-east. This pattern is related to precipitation, where Wyoming rangelands receiving greater than 282 mm of summer precipitation are nearly exclusively grassland, and those below are likely to be shrubland (Driese et al., 1997). The mixed grass prairie in the east (Driese et al., 1997) consists of western wheatgrass (Pascopyrum smithii), blue grama (Bouteloua gracilis), side oats grama (Bouteloua curtipendula), crested wheatgrass (Agropyron cristatum), with locally abundant distributions of cheatgrass (Bromus tectorum) and especially Japanese brome (Bromus japonicus). The High Plains ecoregion is more grass-dominated than the Northwestern Great Plains, with only isolated patches of shrubs. The Middle Rockies, Southern Rockies and Wasatch and Uinta Mountains ecoregions are mostly tree-dominated (Driese et al., 1997) in subalpine areas, with meadows and grasslands dominating alpine areas. These ecoregions also contain patches of mountain sagebrush (Artemisia tridentata ssp. vaseyana). Vegetation cover is typically greater in the mountain ecoregions than surrounding rangelands due to greater precipitation, but growing seasons are much shorter due to cooler temperatures. Overall, the rangelands in our study lie at the ecotone of shrubland and grassland provinces, which have unique climate patterns, disturbance histories, and vegetation communities, and plant architecture.
Fig. 1.
Study area overview. Imagery background is a composite (near-infrared [NIR], red, and blue combination) of Landsat spectral bands concurrent with the 75th percentile of Normalized Difference Vegetation Index (NDVI) values in July 1984–2018. Pixels used as training for Ecological Potential (EP) models are shown in red. Inset map shows the location of the study area in the western U.S. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Methods
Overview
To produce EP component cover maps, we first used Google Earth Engine to produce a composite of the Landsat spectral reflectance associated with the best growing conditions in the archive (1984–2018). The best conditions of each month of the growing season (May-September) were captured separately through obtaining the 75th percentile of NDVI and concurrent Landsat spectral bands throughout the archive (1984–2018). Obtaining data from each growing season month approximates plant phenology, which aids in discrimination among components. To supplement the monthly data, we averaged the Landsat data band-by-band across the months of the growing season to represent total biomass productivity (Gu and Wylie 2010). POLARIS soils data (Chaney et al., 2016) and topographic variables were also selected for inclusion in the component models. We identified ecologically intact sites and extracted data from component cover maps previously generated for ∼2015 conditions (Xian et al., 2015, Rigge et al., 2020a) as training for the EP component cover models. Cubist (RuleQuest Research, 2008) RT models were prepared for each component to generate the EP component cover maps. This method produces EP component cover maps that correspond with the best growing conditions in the least disturbed time period of the Landsat archive. Since the models include 2015 data they do not solely represent potential conditions, but rather the best proxy of these conditions readily available (Rigge et al., 2019). EP component cover maps were not designed to represent 2015, or any particular year, rather the cover expected under the best growing conditions, and in least disturbed condition, in the 1984–2018 Landsat archive. The EP reference then comes dynamically from across space and time on a pixel-by-pixel basis. We regressed the 2015 component cover against the EP component cover and developed maps based on the residuals. Finally, we converted EP and 2015 component cover maps to thematic land cover maps following the method of Rigge et al., (2017). We built on the work of Rigge et al. (2019) by predicting component cover directly, using multiple monthly composites of Landsat imagery across the growing season, generalizing methods to be applicable to all rangelands, and using Landsat (30 m) instead of MODIS (250 m) in model development to address major limitations of spatial resolution. See Supplemental Figure 1 for an overview of processing steps.
Data
Landsat imagery
We leveraged the computational capacity of Google Earth Engine (Gorelick et al., 2017) to identify the 75th percentile values of NDVI from across the entire Landsat archive (1984–2018) in each of the growing season months of May, June, July, August, and September, as a proxy of the best and least disturbed conditions in the study period. The 75th percentile is calculated on a pixel-by-pixel basis, and month-by-month basis, pixel values need not come from the same year. We downloaded the corresponding spectral bands (blue, green, red, near-infrared (NIR), short wave infrared (SWIR) 1, and SWIR 2) in addition to the NDVI (total layers = 7), for each growing season month. Landsat data were Collection 1, Tier 1 surface reflectance from Landsat 5, 7, and 8, and the Quality Assurance (QA) band was used to exclude non-clear (ice, snow, cloud, and shadow) observations in the 75th percentile calculations. The 75th percentile was found to be a good balance between excluding outliers, cloud contamination, and noise found in higher percentile ranks; and still represent ecologically intact conditions. We tested the 90th percentile as used by Rigge et al. (2019) in the Great Basin, which introduced numerous artifacts, inferring that the percentile used may need to be tailored regionally. Final imagery products consisted of (1) 75th percentile of NDVI and the contemporaneous spectral bands in each of five growing season months, and (2) average of each layer value across months in the growing season.
Other independent variables
We selected several key POLARIS soil datasets (Chaney et al., 2016), including available water capacity (m3/m3), clay (%), organic matter (log10%), sand (%), and silt (%) for model inclusions. These data were produced for 6 soil depth increments: 0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm, and 100–200 cm. Next, a 30-m resolution Digital Elevation Model (DEM) was used to calculate aspect, position index, and slope (%). These auxiliary independent variables are critical to improving EP model performance (Herrick et al., 2019, Zhu et al., 2016).
Model development
Training area
To ensure the EP component cover maps represent undisturbed conditions, we trained the EP models on areas of ecologically intact rangeland. First, we excluded pixels identified as non-rangeland (agriculture, urban, water, bodies, and forest; based on spectral index thresholds and U.S. Department of Agriculture, National Agricultural Statistics Service 2013 Cropland Data Layer and National Land Cover Database 2011 agricultural classes) in the 2015 component cover data (Rigge et al., 2020a), within 60 m of roads (TIGER/Line U.S. Census Bureau), and those that burned between 1984–2017 as identified by the Monitoring Trends and Burn Severity (MTBS). Second, since most annual herbaceous cover is non-native, we excluded pixels where annual herbaceous cover composed greater than 20% of the total herbaceous cover. We then placed 25,000 random points within the area meeting the criteria above, extracted pixel values at those random points, and regressed the growing season average NDVI against 2015 bare ground. We developed the exponential regression (bare ground (%) = 295,392 e–0.0072 * growing season average NDVI) to fit 2015 bare ground cover against growing season average NDVI, and convert the former to expected bare ground across the study area (R2 = 0.76, p < 0.05). Since higher bare ground relative to EP is often indicative of degradation, we removed pixels with bare ground cover greater than expected in the exponential regression from the sample pool. Following these procedures, 26.8% of the study area was available to train models (Fig. 1), from which we selected 120,000 random points. Finally, we thinned the number of observations in the mid-range (mean +/– 1 standard deviation) by a factor of 2, which diminishes the tendency to overpredict and underpredict the low and high end of component cover ranges, respectively (Wylie et al., 2019). Those pixels included in the training pool had significantly less annual herbaceous and bare ground, and significantly greater herbaceous, litter, sage, and shrub cover than those excluded from the pool (Table 1). This result suggests that the procedures used to select training data were successful in identifying ecologically intact sites. The approach described above dramatically expands the sample size and geography of potential training area relative to only selecting training data from the high-resolution sites used by Xian et al. (2015) and Rigge et al. (2020a) to produce the 2015 component cover maps. However, error present in the 2015 component cover maps (see Rigge et al., 2020a; Table 2) may be propagated into our training data.
Cubist modeling
We developed Cubist RT (ensemble, boosted approach) models for five dependent variables of EP component cover. We obtained the training data for four of these variables bare ground, litter, sagebrush, and shrub cover from the 2015 component cover products (Xian et al., 2015, Rigge et al., 2020a), available at www.mrlc.gov. For the fifth variable, we chose to predict perennial, as opposed to total herbaceous cover, produced by subtracting annual herbaceous cover from total herbaceous cover. We did not produce EP maps for annual herbaceous cover as it was assumed that most of its constituent species are introduced and therefore EP would be 0%. We extracted the values from both sets of Landsat imagery products (i.e., monthly and growing season average of 75th percentile NDVI and concurrent spectral bands), POLARIS soils datasets, and topographical (aspect, position index, slope) variables to the set of ∼110,000 random points within the training area defined above to serve as independent variables. To parameterize each Cubist model, we used a committee of 5 members (i.e., ensemble) with a maximum of 500 rules and 5% extrapolation allowed. The rules constitute a defined series of multiple piece-wise regression models, and the extrapolation limit allowed predicted values to be 5% beyond the range of training data values. The spatial extent of the EP models was similar to that used in the regional predictions of 2015 components. We calculated the importance of independent variables in each model as the average of model stratification and usage as regression independent variables, which were derived from the model development database. Model stratification reflects how useful each independent variable was in partitioning data space into similar dependent variable responses. The importance value reflects the percentage of pixels in the modeled output influenced by each independent variable.
Departure
From the same set of random points used to train models, we plotted EP component cover against 2015 component cover (i.e., cross validation). A linear regression model was fit for each relationship, all yielding a robust correlation and slope near 1:1 (Table 2). Departures were derived from distance from the overall regression line (regression accounted for minor model biases), as opposed to EP component cover minus 2015 cover. This approach results in the mean regression deviation across all training points to equal ∼0 for each component, at these sites we assume component cover is at or near EP. However, at pixels not included in the training pool, the mean regression deviation (i.e., EP departure) need not have a mean of 0. The resulting departure values are positive in cases where 2015 components are greater than EP component predictions, and negative when EP estimates are greater than 2015 estimates. The departure for each component reflects the difference between EP cover and the 2015 conditions in 2015.
Table 1
Average 2015 component cover in pixels included in training EP component cover maps and those excluded from the training pool. Letters indicate significant differences (p < 0.05) in the mean 2015 component cover between pixels excluded and included in EP training pool, by component (i.e. by column).
Table 2
Linear regression statistics for Ecosystem Potential (EP) component cover regressed against 2015 component cover at training data sites. Departure was derived from the distance from the regression model. All slopes and models are highly significant (p < 0.05).
Land cover
We converted our EP and 2015 component maps into National Land Cover Database cover types using the methods of Rigge et al. (2017). Briefly, this procedure used a series of indices and conditional statements to convert the fractional cover of shrub, bare ground, herbaceous, and litter to barren, shrubland, and grassland land covers. The procedure was slightly modified to accommodate the EP components, specifically we do not consider shrub height since that component has relatively lower accuracy among 2015 components. The departure between EP land cover and 2015 land cover was calculated and summarized by Level III ecoregion. Next, we calculated the proportion of each land cover scenario attributable to known drivers of change, consisting of 1984–2016 MTBS fires, National Wetland Inventory (NWI) wetlands, and Bureau of Land Management vegetation treatments. We included all vegetation treatment types indicated as implemented in the Land Treatment Digital Library compiled by the USGS (Pilliod, 2009). Wetlands were included in the known change as they are subject to change in ground cover (open water, dry [barren] bottom, herbaceous cover, etc.) that are related to water-level fluctuations. Additionally, we calculated the mean 2015 and long-term (1981–2015) Water Year Precipitation (WYPRCP), defined as October-September based on Daymet weather data (Thornton et al., 2016), to serve as strata for analysis.
Statistical summary
We generated EP departure histograms and means by component for the overall study area and all major Level III ecoregions separately. Next, we calculated the mean EP component cover by Level III ecoregion, excluding areas of known change to better isolate weather impacts. Finally, excluding areas of known change, we evaluated the correlations between EP departure, EP component, and 2015 component values with 1981–2015 average WYPRCP, water year minimum temperature (WYTMIN), and water year maximum temperature (WYTMAX).
Weather variables assessed were WYPRCP, WYTMAX, and WYTMIN. We calculated the 1981–2015 average in addition to the yearly averages of 2014 and 2015 for Daymet weather data (Thornton et al., 2016). The years 2014 and 2015 were both found to be considerably wetter and warmer than the long-term average. Therefore, we decided to conduct an analysis investigating the impact of 2014 and 2015 weather conditions on component departures. First, we subtracted long-term weather variables from 2014 and 2015 variables to generate weather departure. Sites with above long-term average WYPRCP in 2015, for example, would be wetter than the long-term average and have a positive departure. Second, we evaluated the correlations between 2014 and 2015 weather departure and component departure. The purpose of this analysis was to determine if EP departure values are more strongly related to 2015 weather conditions or long-term climate. If the latter is the case, then we can be more confident the departure is not ephemeral in nature and merely related to 2014 and 2015 weather.
Validation
Independent validation data were collected in the field in 2015, concurrently with training data for our 2015 component maps (Xian et al., 2015, Rigge et al., 2020a). At each validation point, the fractional cover of all components was observed in a 1 m2 quadrat using ocular estimation every 5 m along a 30 m transect, and the seven observations were averaged. Validation sites (n = 150) were widely distributed across the study area (Fig. 1). At these locations, we calculated the mean component cover from field-measured data, EP component, 2015 component, and EP departure. Next, we calculated the correlation coefficient of EP component and 2015 component regressed on field-measured data. Finally, we conducted an informal visual interpretation of the EP component and departure maps with attention to areas of known disturbances such as MTBS fires, oil and gas development, road and pipeline development, and surface mining clearly visible on circa 2015 Landsat imagery or Google Earth. Our goal in the visual interpretation was to confirm the presence of these features on the EP departure, but not EP component, maps.
Ecological potential limitations and assumptions
Grazing by livestock and wild ungulates represents a potential confounding factor in the selection of training data sites as it is widespread on public (Beschta et al., 2013) and private lands. Impacts of grazing vary based on ecosystem, with some systems relatively tolerant of grazing, however, bare ground generally increases with grazing intensity (Augustine et al., 2012; Porensky et al., 2017). Therefore, we often excluded heavily grazed sites from the training data pool that would be identified as excessive bare ground in the exponential model. Similarly, some light to ungrazed sites may not be sampled, as they could be associated with weedy and invasive species (Porensky et al., 2017). Overall, we propose that our training data represent a light-moderate level of grazing as the basis of comparison.
Another potentially confounding factor is the method in which we prepared imagery to be included in EP models, specifically that we extracted the 75th percentile of NDVI separately for each month. It is unlikely that any given year would have the 75th percentile NDVI across all growing season months. However, the benefit of deriving the 75th percentile NDVI from each month separately, as done here, is to better capture the impacts of phenology. For example, a region could have received plentiful snowfall, resulting in a high May NDVI. This region could also have in the same year (or more likely in a different year), a significant rainfall event in June, resulting in a high July NDVI. These two events may never have occurred in concert in the same year in our temporal record of 1984–2018, but their co-occurrence is within the realm of possibility, especially over a longer time scale. The alternative to this approach would be to look across the entire growing season, and obtain the single highest NDVI value, with no regard given to when in the growing season it occurred, and in the example above miss the double peak, which could influence EP cover estimates. Our premise is that capturing the 75th percentile NDVI of each month separately allows us to better disentangle the timing and duration of NDVI (i.e., greenness) influence on relative EP cover. Moreover, we have found that providing Cubist models more information from across the growing season has improved performance in our base mapping, allowing us to better discriminate among components (Xian et al., 2015, Rigge et al., 2020a). The net result is that our EP cover values for vegetation components are likely somewhat inflated (and conversely bare ground likely deflated). Our goal in EP mapping was to make somewhat optimistic maps of component cover, considering the longest time scales possible with the data available. Our method of calculating EP departure based on regression residuals, as compared to raw difference, removes some of this bias in EP component cover.
Table 3
Average component departure and Ecological Potential (EP) cover by level III ecoregion in unburned areas. Burned area average is based on pixels within 1984–2018 burns.
Results
Spatial patterns
EP components had a higher mean perennial herbaceous, litter, sagebrush, and shrub cover, and lower bare ground than 2015 components (Fig. 2). However, gross spatial patterns of 2015 and EP component cover were similar (Fig. 2). General patterns of EP departure were related to fire history, energy development, surface mining, and vegetation treatments. Burned areas tended to have positive bare ground departure and negative shrub and sagebrush departure since EP components do not capture fire influence, but 2015 components do. Areas of energy development, surface mining, and vegetation treatments were also associated with departure, primarily lower than EP shrub and above EP bare (Supplemental Figure 2). However, the primary factor driving EP departure was weather departure from optimal growing conditions. The EP component cover values were consistent with more productive growing conditions than mapped in 2015 (Table 3). Bare ground cover was higher than EP in the majority (78.8%) of the study area and below EP in 8.2% (Figs. 2,3, Table 3). Conversely, perennial herbaceous, shrub, and sagebrush cover was below EP in 54.4, 58.30, and 53.71% of cases, and above EP in only 29.6, 24.2, and 20.1% of cases, respectively. Removing areas of known change resulted in little difference in ecoregion EP departure means (Table 3), underscoring both the ubiquity of departure and the primary role of weather as a driver of departure. Overall, 2015 bare ground cover was 123.5% of EP (i.e., 23.5% above EP), perennial herbaceous, 88.3%, litter 93.2%, shrub 86.7%, and sagebrush 86.8% of EP.
The EP component predictions successfully removed the influence of burns on EP since the 75th percentile Landsat imagery only captured the influence of burns that occurred before 1984. Within these older burns, EP represented the most productive conditions of the Landsat archive, likely near the end of the 1984–2018 time period. Indeed, visual inspection of the EP component maps showed no influence of burns. The mean EP departure for areas burned between 1984 and 2015 was positive for perennial herbaceous (8.1%), annual herbaceous (3.9%), and bare ground (12.4%), and negative for litter (–1.1%), sagebrush (–2.8%), and shrub (–7.5%) (Table 3). This pattern emphasizes the ability of the EP model to ignore 2015 burns and suggests that burns tend to occur in the most productive portions of the study area. While not analyzed here, Rigge et al. (2019) found that EP departure varied by time since burn in the Great Basin and hence we expect similar patterns in Wyoming.
Climate relationships
Relationships between climate and EP and 2015 component cover were quite similar (Table 4). Although climate data were not inputs of the EP component models, the EP components tended to have a stronger climate relationship than the 2015 components. EP perennial herbaceous, litter and bare ground were most strongly related to WYPRCP, and EP shrub and sagebrush cover were most strongly linked with WYTMAX. Shrub and sagebrush were also strongly negatively related to WYTMIN. Perhaps the most surprising finding was the lack of EP sagebrush relationship with WYPRCP, whereas shrub was positively related to WYPRCP. This finding chiefly resulted from higher WYPRCP and lower sagebrush cover in the east (Fig. 2, Table 3).
Land cover
The land cover departure map (Fig. 4) revealed interesting patterns of potential land cover change. Generally, departure was more widespread in the east, especially in the ecotone between shrubland and grassland. Departures were also common in salt desert scrub basins in the southwest portion of the study area. Localized departures associated with fires, roads, vegetation treatments, mining (especially in the Powder River Basin) and oil and gas development were evident. Overall, 61% of EP barren land cover was barren in 2015, 23% was shrubland, and 16% was grassland (Fig. 5). In 2015, 83% of EP shrubland was present, and 2 and 12% was barren and grassland, respectively. For EP grasslands, 80%, <1%, and 19% were grassland, barren, and shrubland in 2015, respectively. The High Plains ecoregion had the lowest proportion of its area with departure from EP land cover at 7%, the Middle Rockies and Northwestern Great Plains had the greatest at 24% each. Land cover departure in the Middle Rockies was primarily due to EP shrubland existing as grassland in 2015. In the Northwestern Great Plains, 34% of EP shrublands were occupied by grasslands and 9% of EP grasslands were shrubland. In all ecoregions, a greater proportion of EP grasslands existed contemporarily than either EP shrubland or barren.
Fig. 2.
Component cover mapped in contemporarily in 2015 (left) (Xian et al., 2015), Ecosystem Potential (EP) component cover (middle), and EP departure (right). The departure values are positive in cases where 2015 component cover are greater than EP component cover, and negative in reverse.
Fig. 3.
Histograms of Ecological potential (EP) departure for the entire study area (overall) and major level III ecoregions by component. Departure values greater than 0 indicate 2015 cover was greater than EP. See Fig. 1 for delineations of level III ecoregions.
The most prevalent land cover departure was EP shrubland occurring as grassland in 2015, impacting 11.7% of the study area (Table 5). Within this scenario, 9.8% burned during the MTBS record, 2.1% were near wetlands, and 0.7% had vegetation treatment, totaling 12.6% of the area explained by spatially discrete known drivers of change. The second most common land cover departure was EP grassland conversion to shrubland in 2015. Under this scenario, only 5.9% of the change could be attributed to known change, though some additional areas of change may be explained by missing treatments in the Land Treatment Digital Library (Table 5). EP shrubland sites occurring as such in 2015 had significantly higher mean 2015 and long-term WYPRCP than those occurring as barren in 2015. EP shrubland sites occurring as grassland in 2015 had, in turn, higher mean 2015 and long-term WYPCP than those with shrubland land cover in 2015. EP grassland sites echoed these same patterns, with WYPRCP increasing from 2015 barren to shrubland to grassland. Barren land cover existed on <1% of the study area and in general its departures were less stratified by climate than shrubland or grassland.
Table 4
Ecological potential (EP) departure, EP component, and 2015 component spatial correlations (r) with long-term average (from a random sample n = 75,000) water year precipitation (WYPRCP), maximum temperature (WYTMAX), and minimum temperature (WYTMIN). Areas of known change were excluded.
Fig. 4.
Departure between Ecological potential (EP) and 2015 land cover. The ∼280 mm contour of growing season precipitation (GSPRCP) was suggested by Driese et al., (1997) as the approximate shrubland:grassland ecotone. Locations averaging more than 280 mm of GSPCP from 1981–2015 are located north and east of the contour line.
Influence of 2014 and 2015 weather
The weather conditions in 2014 and 2015 were significantly wetter and warmer than the long-term average (Fig. 6). The weather conditions in 2015 and especially 2014 would be expected to influence the 2015 component predictions, and as a result, be the causality of some EP departure. Although these patterns were expected, we wanted to evaluate the degree to which 2014 and 2015 weather departures from the long-term means drove EP departure. Perennial herbaceous EP departure had a weak positive relationship with 2014 and 2015 WYPCP departure but was more strongly related to WYTMIN and WYTMAX departure. Annual herbaceous EP departure had a strong negative correlation with temperature departure in both 2014 and 2015. Shrub and sagebrush EP departure responded similarly to weather departures, both positively related to WYPRCP anomalies and negatively associated with WYTMIN and WYTMAX anomalies (Fig. 7). In general, 2015 weather departures generally had a more significant influence than 2014 departures.
In the portion of the study area with 2015 WRPCP within 25 mm of the long-term mean (4.7% of the study area), bare ground EP departure averaged 7.8%, perennial herbaceous -1.8%, shrub -2.4%, and sagebrush -1.9%. This result, relative to the overall average EP departure in Table 3, indicates that although short-term weather variation drives a portion of departure, the underlying factors are the primary drivers, chief among those being long-term weather variation. Our method of using the residuals to define departure from the EP to 2015 component cover detrends some regional yearly weather influence. For example, mean perennial herbaceous EP departure in principle should not vary if the study area mean WYPRCP was uniformly 100 mm below the long-term mean, as compared to uniformly 100 mm above the mean. However, local anomalies in weather would not be accounted for in the overall regression, and therefore be depicted as EP departure. Although 2014 and 2015 were both wetter than average years, component predictions from 2015 still showed less vegetation cover than EP.
Fig. 5.
Proportional area of 2015 land cover by of Ecological potential (EP) land cover. For example, 82% of the overall area with shrubland EP was shrubland contemporarily, 16% was grassland, and 2% was barren. The horizontal axis is the EP land cover, and color represents 2015 land cover.
Table 5
Land cover departure from EP statistics and drivers. The proportion of study area occupied by each land cover departure scenario, proportion impacted by known drivers of change, and mean 2015 and 1981–2015 Water Year Precipitation (WYPRCP) by scenario. The area column is the proportion of the total study occupied by each scenario. The letters for WYPRCP indicate significant differences (p < 0.05) within a time-period by EP land cover.
Predictor importance
Accuracy assessment of EP component models showed a robust correlation with 2015 components and slopes near 1 (Table 2). Overall NDVI averaged over growing season months was the most important independent variable for EP component development (Supplemental Table 1). For EP litter, August NDVI was most important, and average SWIR band 1 was most important for shrub cover. Although the growing season average Landsat data were most important, the monthly 75th percentile Landsat composites were also important. Monthly spectral data from May were most important to sagebrush, July for shrub cover, August for bare ground and litter, and September for perennial herbaceous (Supplemental Table 2). The interaction of spectral data across time within the growing season (i.e., phenology) is critical for accurate mapping and differentiation of components. Soil depth increments of 0–5 cm and 5–15 cm for all soil properties were most important across EP components, with relative importance decreasing with depth. Clay content was found to be the most crucial soil variable across components, with sand and silt also being quite important. Surprisingly, organic matter (OM) and available water capacity (AWC) were relatively unimportant, possibly resulting from multicollinearity with other variables.
Independent validation
Field observations of component cover were moderately related (r of 0.28 to 0.88) to EP and 2015 component cover (Table 6). EP and 2015 components had similar relationships to field data, but EP relationships were consistently weaker. Mean field-measured values of perennial herbaceous, annual herbaceous, bare ground were higher than both EP and 2015 components, although litter, sagebrush, and shrub were lower. This finding suggests that there was some bias in the 2015 components relative to field observations, resulting in bias in the EP components. Moreover, this finding indicates that some cases of departure could be caused by errors in both the 2015 and EP component maps. However, EP component values reflected more productive conditions than 2015 components. EP components had a higher mean perennial herbaceous, litter, sagebrush, and shrub cover, and lower bare ground relative to field-measured component cover (Table 6). Similar EP modeling undertaken by Rigge et al. (2019) in the Great Basin was successful in documenting the widespread and lasting impacts of fire on vegetation structure. Moreover, Table 3 demonstrates that burned areas on average had more positive perennial herbaceous, annual herbaceous and bare ground, and more negative sagebrush and shrub EP departures than unburned areas, though Porensky and Blumenthal (2016) found that fire may not promote cheatgrass invasion in the western Great Plains.
Fig. 6.
Study area mean and standard deviation for the 2014, 2015, and 1981–2015 average of A) Water Year Precipitation (WYPRCP), B) Water Year Maximum Temperature (WYTMAX), and C) Water Year Minimum Temperature (WYTMIN). Letters indicate significant differences (p < 0.05) among the time period mean climate variable, within each panel (i.e. within each climate variable).
Fig. 7.
Spatial relationships (r) between EP departure and 2014 and 2015 weather departure from 1981–2015 mean by component (from a random sample n = 75,000). Results are shown for water year precipitation (WYPRCP), maximum temperature (WYTMAX), and minimum temperature (WYTMIN). All correlations are significant (p < 0.05).
Table 6
Mean component cover at independent validation sites (n = 150) for field measured data, Ecological potential (EP) departure, EP component, and 2015 component cover. Correlation (r) of field measured data against EP component and 2015 component cover. Asterisks indicate significant (p < 0.05) correlations.
Discussion
The approach used to map EP and departure documented in the current work improved upon our previous work in several critical areas. First, Rigge et al. (2019) used MODIS growing season averaged NDVI as the dependent variable of EP modeling, which does not discriminate among vegetation types as effectively as monthly data. In addition to NDVI, we included spectral values from multiple bands for all growing season months, and these additional inputs were valuable in EP model development. Moreover, Reeves and Bagget (2014) found that the disturbance grain size of riparian degradation and other anthropogenic influence was too small to be resolved with MODIS data; therefore, mapping at the Landsat resolution (30 m) represents a critical improvement. Second, we use growing season average NDVI as only one of several independent variables to predict EP component cover using RT models, whereas Rigge et al. (2019) used regression functions to convert growing season average NDVI to EP component cover directly. Third, predicting components using RT models allowed us to avoid the usage of component cover weighting factors employed by Rigge et al. (2019), which rated the somewhat subjective value of different components in the composition of a healthy sagebrush steppe. These process improvements yielded more robust cross validation (Table 2) relative to the previous work. Independent validation data (Table 6) show that EP components were reflective of more productive conditions than experienced in 2015, a wet year. The weaker relationships with field-measured data in EP components as compared to 2015 components was expected, since disturbance, management, and less productive weather conditions in 2015 relative to the EP would cause 2015 values to deviate from EP.
EP component models were mostly driven by spectral data, especially the growing season average NDVI (Supplemental Table 1). In general, growing season-averaged Landsat data were more important than monthly data, but of the monthly data, August and September were most significant (Supplemental Table 2). Late summer tends to be the driest portion of the growing season in the study area, leading to greater landscape stratification, and improved spectral discrimination (Vande Kamp et al., 2013). The role of non-spectral data was lower than spectral data, although both EP and 2015 component predictions were strongly correlated with soil properties, implying that much of the relevant information to mapping EP components was found in spectral data alone. Clay content was most important among soil variables, and shallow depths most important overall. This was somewhat surprising for shrubs and sagebrush since sagebrush can access deeper moisture pools, at least 91 cm (Sturges 1977) but can be deeper. However, Sturges (1977) also found that maximum lateral root spread of sagebrush at a Wyoming site was in the surficial 30.5 cm, especially in the top 15 cm, perhaps reflecting the importance of these depths in our EP component models. Moreover, deeper roots were found negligible to overall plant water dynamics
Spatial patterns
Our EP component and departure maps successfully capture management-relevant spatial patterns of both within-state change and land-cover change. The EP component models predicted reasonable patterns (Fig. 2), consistent with ecological expectations. For example, EP shrub and sagebrush cover is higher in the west and EP bare ground is negatively associated with WYPRCP. From a regional perspective, positive bare ground departure was most pervasive and significant, resulting from negative departure in perennial herbaceous, litter, and shrub cover (Fig. 2, Table 3). Spatially discrete disturbances of burns, energy development, and vegetation treatment were clearly distinguishable on EP departure maps through visual interpretation (Supplemental Figure 2). Reeves and Bagget (2014) found that oil and gas development, livestock watering facilities, and prairie dog (Cynomys spp.) towns were primary degradation agents in the Great Plains. Allred et al. (2015) reported that an average of 50,000 new oil wells was constructed each year in central North America since 2000, estimated as the largest driver of land use change in the U.S. (Trainor et al., 2016). In Wyoming, energy development has been ongoing for over 100 years, but a substantial portion has occurred from circa 2000 onwards (Allred et al., 2015), which will appear as EP departures. This development has long-term impacts on vegetation structure. For example, Avirmed et al. (2015) estimated that 87 years are required for a full recovery of big sagebrush cover in abandoned oil pads, and forbs may never recover. Much of the vegetation that does occur on abandoned oil pads are weedy annuals (Waller et al., 2018) or grasses and non-sagebrush shrubs (Avirmed et al., 2015). Although we found considerable amounts of disturbed lands, Wyoming was found to have a relatively low proportion of degradation compared to surrounding areas (Herrick et al., 2010).
The predominant pattern in locations with no known disturbance was below EP cover for vegetation components and above EP for bare ground. This pattern could arise due to a combination of grazing, management, weather impacts, and potential bias in the EP models. While our training data set does include sites with light-moderate levels of grazing, sites with heavy grazing would commonly have higher cover of bare ground (Augustine et al., 2012, Porensky et al., 2017) and therefore be excluded from the training data pool. Thus sites with moderate-heavy grazing would commonly induce positive departures in EP bare ground and negative departure in EP vegetation components. Weather in 2015 was wetter and warmer than the study period average (Fig. 6). However, other years in the Landsat record were much wetter. Average WYPRCP in 1995, the wettest year in the time-series, was 500.1 mm +/– 112 mm, considerably wetter than in 2015. Moreover, a severe drought occurred in 2012, with mean WYPCP of 233.8 mm +/– 87.2. The severe drought of 2012 and long-term drought of 2000–2010 may also have lagged effects (Moran et al., 2014), resulting in lower 2015 vegetation cover and therefore EP departure. The drier conditions in 2015 relative to the maximum WYPRCP observed in the climate record and disturbance are reflected in the direction and magnitude of EP departure (Figs. 2 and 3).
EP and 2015 component climate relationships
EP and 2015 component cover spatial relationships with climate are generally as expected (Table 4) with positive associations between WYPRCP and vegetation components, and differing relationship by component. EP departure associations were somewhat more difficult to interpret, as sagebrush, for example, has a diverse response to weather patterns (Monroe et al., 2020) and disturbance (e.g., silver sagebrush [Artemisia cana] (White and Currie, 1983) compared to big sagebrush (Baker, 2006)). Unlike results from the Great Basin (Shi et al., 2018), we found both EP and 2015 sagebrush cover was negatively related to WYTMIN, and the sagebrush EP departure was positively associated with warmer WTYMIN and WYTMAX (Table 4). This could result from differences in the stratification of vegetation types by elevation in the Great Basin as compared to Wyoming, where higher elevation (i.e., cooler TMIN) sites in the Great Basin tend to be drier and have more sagebrush cover than those in Wyoming, which are more forest dominant. Moreover, shrub and sagebrush EP and 2015 cover were more strongly related to WTMAX and WYTMIN than to WYPRCP. Our study area is on average wetter and cooler than the Great Basin. Thus, sagebrush may be more cold temperature limited as opposed to water limited in the study area. Some error in the climate relationships could be related to known tendencies to overpredict and underpredict component cover at the low and high ends of the histogram, respectively. This pattern would be consistent between the EP and 2015 components, however, suggesting minimal impacts to EP departure values.
Land cover
Relative to an early Landsat-based map of Wyoming vegetation types associated with gap analysis (Driese et al., 1997), we show more EP shrubland in the Laramie Valley and east slopes of the Big Horn Mountains and more grassland in east-central Wyoming. Because of disturbances resetting land cover to any earlier successional state, the primary direction of land cover departure was EP shrubland occurring as grassland in 2015 (Figs. 4 and 5). On the other hand, EP grassland occurrence as shrubland in 2015 was the second most common departure. Similar to the woody plant expansion of chiefly Juniperus spp. into mesic grasslands (Briggs et al., 2005), the expansion of the shrubland system into EP grasslands may be due to reduced fire frequency in the Northern Great Plains relative to the historical 5- to 10-years return interval at most locations (Sieg, 1997, Wendtland and Dodd, 1992, Wright and Bailey, 1982,). Conversely, when looking at only large fires, Dennison et al. (2014) reported an increasing trend in the Northern Great Plains. In the absence of fire, sagebrush cover increases on uplands in the western portion of the Northern Great Plains (i.e., near the ecotone between the Northwestern Great Plains/High Plains and the Wyoming Basin ecorgions) (Wright and Bailey, 1982). The primary concentration of land cover departure is near the ecotone between shrubland and grassland in the eastern third of Wyoming. Like the rest of Wyoming, MTBS data reveal limited recent fire activity near the ecotone and known drivers of change explained little EP land cover departure (Table 5)
Climate change can also drive movement in ecotones, which are especially sensitive to climatic changes (Peters, 2002). Driese et al. (1997) determined the ∼280-mm contour of average growing season precipitation as the approximate boundary of the shrubland:grassland ecotone in Wyoming. This contour generally lies to the east of the region of greatest EP land cover departure (Fig. 4). The ecotone is very sensitive to thresholds in the Rigge et al. (2017) model and the generally wetter than average conditions in 2015 would be expected to push the ecotone west as we observed. Moreover, local edaphic conditions were found to be more important than regional climate trends in resolving the ecotone (Driese et al., 1997). Overall, our data suggest the ecotone is more diffuse in 2015 than EP, indicative of management as a driver. We must also note that the complex mixture of grasses and shrubs near the ecotone often produce a complex spectral environment which can increase mapping error in both the EP and 2015 components (Rigge et al., 2020a).
Ecological potential
Bai et al. (2008) defined land degradation as long-term decline in ecosystem function and productivity. Degraded and early successional states are frequently stable, and common in the sagebrush steppe (Laycock, 1991). Degradation often results in reduced and/or altered vegetation cover and loss of soil organic carbon through erosion (Van der Esch et al., 2017). Our EP departure maps provide a quantitative resource in evaluating spatial patterns of departure from potential conditions observed over the Landsat record. The lack of a strong spatial relationship between 2015 weather departures and EP departure (Fig. 7) reveals that our methods partly control for the effects of climate. In other words, the EP departure results are somewhat stabilized from year-to-year weather variation, and instead are more strongly related to long-term perturbations. This property of dampened sensitivity to short-term weather variation is beneficial in discriminating the long-term impacts of land management (Van der Esch et al., 2017). In undisturbed areas, the EP component cover represents the conditions of the 2015 plant community under ideal conditions. In areas disturbed either during or before the Landsat record, the EP components predict the component cover concomitant with the best growing conditions from the portion of the archive with the least disturbance; this can be both prior to the disturbance occurring or reflective of the most complete recovery. Thus, EP components in disturbed areas do not represent the EP precisely, but the nearest to that which was observed in the Landsat archive.
Implications
We use time as the potential condition, dynamically selected through space, from the Landsat record to correspond with the least disturbed and best growing conditions. By predicting EP component cover separately, we avoid the uniform treatment of the study area in regard to disentailing component cover from NDVI. Combining EP fractional components could yield an impression of potential, and EP departure combinations among components could suggest current state. For example, Davies et al. (2012) proposed a state-and-transition model based on two axes, 1) shrub-herbaceous ratio and 2) native-invasive herbaceous ratio. The shrub-herbaceous axis is important at a community scale as it serves as a surrogate for the recolonization likelihood/potential of shrubs. The native-invasive axis describes the character of the herbaceous presence and thereby the likelihood, role, and impact of fire. The combination of the two axes describes the 1) site intactness, 2) likelihood of recovery, 3) resistance and resilience. Our EP and 2015 component cover data could be employed in such a method. Moreover, both the EP and 2015 component data could be split into states on an existing Natural Resource Conservation Service state-and-transition model. One key problem with this method is that the exact thresholds between states (tipping points) are difficult to define as they represent complex interactions among climate, disturbance, and existing vegetation cover (Briske et al., 2005). Moreover, functional thresholds will lag structural thresholds (Briske et al., 2005) defined by EP departure. Another issue is that fires in grass-dominated communities often do not result in state change, but they do prevent colonization by shrubs (Davies et al., 2012), complicating state-and-transition modeling. A simpler alternative may be monitoring change in patch size or vegetation biomass, represented by perennial herbaceous and shrub EP departure, which can indicate that state change could soon occur (Bestelmeyer et al., 2011). Another option is to examine component EP departure in concert, where a particular state change may be present as a series of departure values (e.g., positive values for bare ground and perennial herbaceous, negative for shrub, etc.). Moreover, pixels with EP departures in multiple components can be interpreted as providing higher confidence in the departure.
Perryman et al. (2018) argue that the goal of returning Great Basin ecosystems to pre-settlement conditions has failed as altered environmental conditions and cheatgrass colonization have radically altered the landscape. Climate change is likely to alter landscape condition further, with an increase in the amount and spatial variance of sagebrush cover expected in southwestern Wyoming (Tredennick et al., 2016). Conversely, Homer et al. (2015) found an increasing trend in bare ground from 1984–2011 and decreasing trend in shrub, herbaceous, litter, and sagebrush cover. We find that the study area landscape has been altered, and net EP departure patterns are largely consistent with the component cover trends reported by Homer et al. (2015). The net decrease in sagebrush cover, considering the climate-driven expectation of increased cover (Tredennick et al., 2016), implies that disturbance was a more significant driver of change to that component over the study period than climate change. As restoration to pre-settlement is often not attainable, our usage of the wettest and least disturbed conditions observed in the Landsat record as a proxy for EP may be a more appropriate baseline than pre-settlement vegetation type and condition.
Conclusions
The approach documented in this paper was designed to apply to any rangeland ecosystem in which component cover data are available. Moreover, selecting spectral data associated with different NDVI percentile values or different time periods in Google Earth Engine will result in unique products with a range of applications. Some EP departure is related to the relative spatial patterns of 2014 and 2015 weather conditions, but absolute difference in 2014 and 2015 weather conditions from the long-term climate reference (i.e., best growing conditions represented by the 75th percentile NDVI across the Landsat archive) is the primary driver. Although growing season average NDVI was the most important independent variable across components, we also document the importance of spectral data from across the growing season, in addition to soils data.
EP and 2015 component cover were strongly related to climate and generally followed the expected regional patterns. The 2015 conditions in the study area showed reduced shrub, sagebrush, litter, and perennial herbaceous cover and greater bare ground relative to EP. Known disturbances, such as energy development, fires, and vegetation treatments were clearly visible on the departure maps. No 2015 disturbances were visible in the EP component maps, supporting the validity of departure maps. The most frequent departure from EP land cover was shrubland conversion to grassland, secondary in importance was grassland conversion to shrubland. The former can be explained only in small part by known disturbance, especially fire. The lack of known drivers implies that these land-cover changes were associated with climate change, management, and in some cases, scarcity of fire. These land-cover departures have broadened the ecotone between shrubland and grassland relative to EP. Future work will investigate patterns related to EP departure in the context of resiliency and explore the usage of EP components in a state-and-transition model. EP data are available in Rigge et al. (2020b).
Acknowledgments
The authors would like to thank the Bureau of Land Management and the USGS for their funding support of this research. Work by Rigge and Shi was performed under USGS contract 140G0119C0001. We acknowledge Devendra Dahal for producing Google Earth Engine composites. We also thank Todd Hawbaker, Michael O'Donnell, and two anonymous reviewers for reviewing drafts of this manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.rama.2020.03.009.
