Projected change in rangeland fractional component cover across the sagebrush biome under climate change through 2085

. Climate change over the past century has altered vegetation community composition and species distributions across rangelands in the western United States. The scale and magnitude of climatic in ﬂ uences are unknown. While many studies have projected the effects of climate change using several modeling approaches, none has evaluated the impacts to fractional component cover at a 30-m resolution across the full sagebrush ( Artemisia spp.) biome. We used fractional component cover data for rangeland functional groups and weather data from the 1985 to 2018 reference period in conjunction with soils and topography data to develop empirical models describing the spatiotemporal variation in component cover. To investigate the rami ﬁ cations of future change across the western United States, we extended models based on historical relationships over the reference period to model landscape effects based on future weather conditions from two emission scenarios and three time periods (2020s, 2050s, and 2080s). We tested both generalized additive models (GAMs) and regression tree models, ﬁ nding that the former led to superior spatial and statistical results. Our results indicate more xeric vegetation across most of the study area, with an increasing dominance of non-sagebrush shrubs, annual herbaceous cover, and bare ground over herbaceous and sagebrush cover in both the representative concentration pathway (RCP) 4.5 and 8.5 scenarios. In general, both scenarios yielded similar results, but RCP 8.5 tended to be more extreme, with greater change relative to the reference period. Results demonstrate that in cool sites some degree of warming to growing season maximum temperature or nongrowing season minimum temperature could be bene ﬁ cial to sagebrush and shrub growth. However, warming nongrowing season maximum temperature was bene ﬁ cial to shrub, but not to sagebrush growth. Our results inform rangeland managers of potential future vegetation composition, cover, and species distributions, which could improve prioritization of conservation and restoration efforts.


INTRODUCTION
Climate change has altered plant communities substantially in many locations globally (Polley et al. 2013). These changes are associated with trends in temperature, precipitation, and atmospheric CO 2 , which can alter plant biochemistry and physiology (Lucht et al. 2006). Climate projections indicate community composition and function (Izaurralde et al. 2011), and species ranges will continue to shift in the future with changes to the amount and timing of precipitation, soil moisture, evapotranspiration, and occurrence of drought (Collins et al. 2013, Polley et al. 2013. Moreover, many climate models indicate increased variability in weather (Debinski et al. 2010, Collins et al. 2013, which is already quite variable in western U.S. rangelands (Reeves et al. 2020). The nature of precipitation events is expected to change, often to a lower number of more intense events (Collins et al. 2013), which is problematic in dryland environments where plants are vulnerable to increased aridity (Munson 2013, Palmquist et al. 2016b. Ecohydrological response to climate change is likely to vary spatially across current climate gradients and vegetation types (Palmquist et al. 2016a). Sagebrush (Artemisia spp.), for example, is a widespread shrub of ecological importance across the American West (Davies et al. 2011, Renwick et al. 2017. The widespread occurrence of sagebrush across range of climatic and biophysical conditions is related to its many adaptive features (Lysne 2005). However, ongoing (Klos et al. 2014) and expected continued reductions to snowpack and earlier melt (Harte et al. 2015) may shift the timing of the hydrological season to early in the year (Barnett et al. 2005, Klos et al. 2014, Palmquist et al. 2016a, often resulting in water stress across more of the growing season (Flerchinger et al. 2019). Projections of overall increased precipitation in the future may not offset the impacts to sagebrush resulting from the loss of water storage capacity in snow and projected warming (Barnett et al. 2005). Although sagebrush individuals may survive extremely dry soil conditions in the short term, leaf area and canopy cover are greatly reduced (Kolb and Sperry 1999). In many cases, these hydrological changes will decrease the likelihood of seedling survival, key to maintaining a healthy population, in the future (Lysne 2005, Bradford et al. 2014, Palmquist et al. 2016b), especially in drier sites (Renwick et al. 2017, Kleinhesselink and. At the same time, recruitment and growth may increase in the cooler (Renwick et al. 2017, Kleinhesselink and and wetter (Flerchinger et al. 2019) portions of the current sagebrush range.
Rangeland managers require information on the condition of natural resources over large areas and time scales and information on potential climate change impacts (Soulard and Sleeter 2012, Creutzberg et al. 2015, Munson et al. 2016a. Managers need to make decisions now that will impact the landscape for decades into the future (Ford et al. 2019). Even without the uncertainties posed by climate change, restoration of dryland ecosystems is challenging (Shriver et al. 2018). Recently, a new generation of modeling procedures has been used to develop yearly fractional component (i.e., 0-100% cover of shrubs, bare ground, herbs) cover maps of rangelands across large expanses of the western United States through the Landsat archive from 1984 to 2018 (Homer et al. 2013, Jones et al. 2018, Shi et al. 2018, Rigge et al. 2019a, 2021. These data reflect annual change in functional group abundance in relation to interannual variation in weather, disturbance, management, and succession (Rigge et al. 2019a;Reeves et al. 2020).
Researchers have used a variety of modeling approaches to investigate the potential impacts of climatic change on rangeland ecosystem structure and health. Still and Richardson (2015) used random forests to model the climate envelope of Wyoming big sagebrush (Artemisia tridentata ssp. wyomingensis) through 2050, predicting a 39% loss in suitable climate. Reeves and Bagne (2016) projected climate change impacts across U.S. rangelands, finding an increase in forage production in the north, overall shift to more herbaceous cover (and less shrub), and more interannual variability in forage production in most areas, especially in Wyoming and Nevada. Creutzberg et al. (2015) used a dynamic global vegetation model, climate envelope model, state-and-transition simulation model (STSM), and a sage grouse habitat model combined with climate change and management scenarios to investigate the potential future condition of southeastern Oregon rangeland. Population models can be used to model the climate change impacts to sagebrush, but due to the nature of such observations and computational demands in modeling, population models tend to be applied at local scales (Tredennick et al. 2016). Another approach is to apply empirical climate relationships observed in historical data with simulation software to evaluate the impact of diverging climate scenarios (Soulard and Rigge 2020). Finally, Homer et al. (2015) extended linear correlations between rangeland fractional component cover and historical precipitation data to project cover given two climate scenarios in 2050.
We use U.S. Geological Survey (USGS) backin-time (BIT) fractional component cover data for rangeland functional groups (sagebrush, shrub, herbaceous, annual herbaceous, litter, and bare ground) and climate data from 1985 to 2018 (Rigge et al. 2019a(Rigge et al. , 2021 in conjunction with soils and topography data to develop models describing the empirical spatiotemporal variation in cover across the sagebrush biome of the western United States. Our primary objective is to apply the models developed over the 1985-2018 reference period to future weather conditions based on two emission scenarios (representative concentration pathway [RCP] 4.5 and RCP 8.5) and three time periods (2020s, 2050s, and 2080s) to evaluate change in functional group abundance within the respective range of each across unprecedented extent. We then model component cover under these scenarios and evaluate spatial and temporal patterns of change relative to the reference period. We test two modeling methods to determine which produces the most accurate results. Because the interaction of disturbance, interannual variation in climate, soils, and interacting effects is needed to fully understand the drivers of vegetation change (Morris et al. 2016), we included soils and topographic data as spatial modifiers of future projections (Munson et al. 2016b). Our 30-m projections of component cover are based on empirical spatiotemporally informed climate-biophysical relationships developed across 33 yr in the entire sagebrush biome. This represents a major improvement upon previous modeling efforts in terms of the robustness and comprehensiveness of training data and spatial resolution of model outputs. Our projections inform rangeland managers of potential future vegetation composition, cover, and species distributions, which could improve prioritization of conservation and restoration efforts.

Study area
We investigated the sagebrush biome of the western United States identified by Jeffries and Finn (2019). A total of 1,760,134 km 2 is included in this region, of which 1,302,499 km 2 was identified as rangeland. Elevation ranged from 66 to 4210 m, averaging 1670 m. Water year mean precipitation across the sagebrush biome ranges from 10 to 2360 mm, with a mean of 324 mm. Mean annual minimum and maximum temperatures are 0.15°and 14.38°C, respectively. Climate data were based on 1985-2018 means from PRISM (Parameter-elevation Relationships on Independent Slopes Model) data (Daly et al. 2007

Overview
For the period of 1985-2018, we developed empirical spatiotemporal relationships between annual rangeland fractional component data, weather data, topography, and soil variables. These relationships were obtained by testing two different models (cubist regression tree [RT], RuleQuest Research 2008;and generalized additive models [GAMs]). Both models were trained in rangeland portions of the sagebrush biome, and across the entire time series, in pixels that were unburned throughout the study period and where no known vegetation treatments occurred (i.e., the reference period). Next, we applied the models developed for the range of conditions observed in the reference period to projected future climate data based on two emission scenarios and three time periods. Because climate data are the only model inputs that change between reference period and future conditions, they are the sole driver of change in the future, though the soils and topography modulate the response.

Dependent variables
The primary dependent variable in our study is the BIT time-series fractional component data from 1985 to 2018 (Rigge et al. 2021). The BIT data were themselves trained based on fractional rangeland component cover maps developed across the western United States representing circa 2016 conditions . This dataset was developed through extensive field observations of component cover collected on high-resolution (2 m) satellite images (n = 331) including WorldView-2, QuickBird, and Pleiades. Cubist regression trees (RTs) were used to predict component cover over high-resolution v www.esajournals.org image extents. Component cover predictions were degraded to 30 m to serve as training for regional predictions, which included three seasons of Landsat imagery, topographic variables, and spectral indices as inputs. Again, cubist RT models were used to predict component cover at a regional scale. Non-rangeland areas including urban, water, agriculture, and forest cover were identified using the National Land Cover Database, spectral index thresholds, and the 2013 Cropland Data Layer from the U.S. Department of Agriculture, National Agricultural Statistics Service, and excluded from mapping products. Products were harmonized by ensuring component relationships were intact (e.g., shrub cover > sagebrush cover), and were mosaicked across regions.
A series of modeling procedures conceived of by Homer et al. (2013) and further developed by Shi et al. (2018) and Rigge et al. (2019aRigge et al. ( , b, 2021 was used to predict annual fractional component cover across the Landsat archive from 1985 to 2018, using the circa 2016 maps as training. This method focuses on changes observed in Landsat imagery and is described in detail by Rigge et al. (2019Rigge et al. ( , 2021. The basic structure of the process involves using change vector and change fraction methods to screen for spectral change on a bandby-band basis between the summer and fall image of each year and the summer and fall image of the year used to produce the circa 2016 base map. Change identified in both summer and fall images was selected as the final changed area in which fractional component cover is likely to differ from the value of the base year. Training data points were randomly located within the area identified as non-changed relative to the base year and pooled together across space and time for each component. We selected several independent variables for inclusion in modeling, and the summer and fall Landsat image pixel values from six spectral bands, spectral indices, and topographic variables. The component cover training data and independent variables were input into cubist RT models, which were developed spatiotemporally but applied to individual years. Pixels flagged as change relative to the base year were labeled with the component cover prediction, and the base year value was maintained elsewhere. A series of post-processing steps removed noise and illogical change in the pixel time series and ensured accurate postburn trajectories. The final BIT dataset consists of annual maps of shrub, sagebrush, herbaceous, annual herbaceous, litter, and bare ground from 1985 to 2018, mosaicked across the sagebrush biome. The BIT data have been robustly validated by Rigge et al. (2019b) who compared BIT data to a series of predictions made on high-resolution satellite imagery finding strong temporal relationships with a mean coefficient of determinations (R 2 ) of 0.63 and root-mean-square error (RMSE) of 5.47% and strong spatiotemporal correlations with a mean R 2 of 0.52 and RMSE of 7.89% across components. Data collected at longterm monitoring plots (n = 126) over a 10-yr period of 2006-2018 in southwestern Wyoming were compared to BIT data, again showing robust temporal and spatiotemporal correlations, with a mean R 2 of 0.46 across components . BIT data are available for download , and a web application to visualize data and download custom extents is found at www.mrlc.gov/data.

Independent variables
We used 30-m gridded POLARIS soil data (Chaney et al. 2016) for available water capacity (AWC), clay, organic matter (OM), sand, and silt for two depth increments: 0-30 cm and 31-200 cm. Topographic variables of aspect, elevation, position index, and slope were derived from a 30-m digital elevation model (DEM). We converted aspect from polar (i.e., bearing) to Cartesian coordinates (e.g., 90°[east] has coordinates of x: 0, y: 1). Topographic variables are critical to define the landscape position of sites which combined with soil texture is intrinsically linked to soil water availability (Debinski et al. 2010, Duniway et al. 2010. Using 800-m gridded PRISM data, we calculated the growing season (GS; April-September) and nongrowing season (NGS; October-March) average minimum (MIN) and maximum (Max) temperature (T), and the GS and NGS precipitation (PRCP) totals for each year from 1984 to 2018. This process yielded six annual climate variables GSPRCP, NGSPRCP, GSTMAX, NGSTMAX, GSTMIN, and NGSTMIN. Our summarization of six variables sufficiently captures the critical dynamics of NGSTMIN and v www.esajournals.org NGSTMAX (Shi et al. 2018). In much of the sagebrush biome, precipitation is concentrated in the winter, making snowpack and overall winter moisture, a major crux for sagebrush propagation (Palmquist et al. 2016b) and limitation to the primary growing season of March-June. In the north and east portions of the biome, precipitation tends to less winter-dominated, where precipitation patterns in April-July are strongly related to net primary production (NPP, e.g., Chen et al. 2019). However, future shifts in precipitation and temperature seasonality vary through space and time (Lucht et al. 2006), and there is a high degree of uncertainty in estimates of summer drought (Creutzberg et al. 2015). We chose to create broad climate integration windows (GS and NGS) to capture seasonality, while creating parsimonious models and readily interpretable results and limiting the influence of future projection uncertainty in seasonality.
For the climate projections, we used monthly data downscaled using bilinear interpolation to 1000 m across the western United States. Projections were from an ensemble of 15 coupled model intercomparison project phase 5 atmosphere-ocean general circulation models (Adapt-West Project 2015). We obtained data for the low emission scenarios RCP 4.5 (stabilized emissions with 4.5 W/m 2 of radiative forcing and high emissions) and RCP 8.5 (increased emissions through early to midcentury with 8.5 W/m 2 of radiative forcing). Data were available for three time periods, 2020s, 2050s, and 2080s, representing the mean of 2011-2040, 2041-2070, and 2071-2100, respectively. Using these data, we constructed the same set of six climate variables used in the reference period. Reference period and projected climate data were resampled to 30m resolution using nearest neighbor interpolation to match that of other input data.

Model development
We tested two models for each component; cubist RT and GAM. A GAM relates a univariate response to a set of independent variables, merging the properties of generalized linear models (GLMs) and additive models (Hastie andTibshirani, 1986, Wood 2017). One fundamental challenge of GAMs is to set the amount of function smoothness, set by the number of degrees of freedom, and the penalty of the wiggliness of the model. How to estimate these parameters and how to practically represent the smooth functions are the main statistical questions introduced by moving from GLM to GAM.
We excluded areas identified by MTBS as burned and those with implemented vegetation treatments identified by the Land Treatment Digital Library (Pilliod 2009) from 1985 to 2018 from the training data pool, a total of 14.6% of the study area. While fire is a major driver of change in the study area, including fire or vegetation treatment influenced sites in the training data pool would corrupt the climate signal. The relationships between climate variation and component cover in the reference period would no longer reflect the sole influence of climate. Moreover, Rigge et al. (2021) found that 92% of the BIT pixels had component change over the 1985-2018 reference period. Given that 14.6% of the study area had a known change driver (fire and vegetation treatments), the driver of change was unknown in~77% of the study area. Our fundamental goal in this paper was to capture the climate relationships observed in these pixels in the reference period and apply them to the future.
Training data were prepared by placing random points (n = 90,000 per year, total n = 2,880,000) within the rangeland portions of the sagebrush biome ( Fig. 1). Point locations were randomized across years. Independent data variables (n = 20) were extracted to these points: (1) 1985-2018 climate data for six variables (GSPRCP, NGSPRCP, GSTMIN, NGSTMIN, GSTMAX, and NGSTMAX), (2) four topographic variables, and (3) five soil variables for two depth increments each. Additionally, we extract the dependent variables, the fractional cover of six components, to the points. Collecting training data across space and time avoids some of the limitations of spatial and temporal models (Renwick et al. 2017).
We have designed our model workflow and training data preparation to avoid issues related to spatial or temporal autocorrelation. Our overarching goal was to create a model for each component to describe the general relationships observed through space and time to independent variables that would be readily transferable to future climate conditions. Specifically, we avoid autocorrelation through randomizing point locations across years and excluding spatial (i.e., latitude and longitude) and temporal (i.e., year) v www.esajournals.org data in model development. The random point locations for each year were attributed with data from only that year. In other words, our models for each component are agnostic to space and time, rather they capitalize on relationships observed across all space and time, without specifically accounting for their effect.

Model selection
We evaluated GAM and RT model based on their ability to minimize residuals and have the most robust relationships with training data. Specifically, we generated GAM and RT predictions for 1990, 2000, and 2010 compared them to the corresponding yearly BIT predictions. For the GAMs, we tested several degrees of freedom options: k = 40, 6, and 3. Among the GAMs, six degrees of freedom produced the strongest models with the least amount of spurious relationships included (observed in the k = 40 models), while not oversimplifying relationships (observed in the k = 3 models), these were selected as the final models. We inserted weather data from 1990, 2000, and 2010, to conduct cross-validation across a range of climate conditions in the reference period (Table 1). Relative to the BIT predictions, we found that GAM predictions were superior to RT in terms of having mean predicted values, standard deviation, and range of values more similar to BIT (Table 1), and based on expert opinion, a better characterization of overall landscape patterns. Model fit of final GAMs relative to BIT training data is presented in Appendix S1: Fig. S1.

Future component cover prediction
Based on the final selected GAM developed for each component in the reference period, we replaced climate data with future projections for the 2020s, 2050s, and 2080s, independently. We then predicted future component cover spatially at 30-m resolution across the study area based on the empirical relationships observed in space and time of the reference period. We evaluated the occurrence and degree of extrapolation to novel climate conditions in future projections relative to the reference period; novel is defined as conditions that were never observed in the reference period throughout the sagebrush biome. This is critical as spatial and temporal relationships between biomass production (or cover) and precipitation determined with ground or remote sensing data should not be extrapolated (Munson et al. 2019). Relative to the climate conditions observed through time and space in the reference period, we found very little novel climate conditions, even in the most extreme climate scenario of RCP 8.5 in the 2080s (Fig. 2). This demonstrates (1) the vast amount variation in space and time in the reference period, and (2) model extrapolation to novel climate conditions is limited.
Following the completion of all GAM predictions, we implemented a post-processing procedure on each set of component predictions for each year/scenario. This process ensured the sum of primary components (bare ground, herbaceous, litter, and shrub cover) equaled 100% by keeping the original modeled proportions intact. Next, we ensured that shrub cover was greater than sagebrush cover, and herbaceous cover was greater than annual herbaceous cover. Projected future component cover layers are available at Rigge (2020).

Model limitations
Our study has several caveats and limitations. First, the distribution of precipitation in terms of number and magnitude of events and type of precipitation (snow vs. rain) is important, but beyond the scope of this paper. Second, the effects of elevated CO 2 itself are largely positive on plant growth, especially in regions in which plant growth is moderately water-limited, and in C 3 plants (Izaurralde et al. 2011, Polley et al. 2013). We do not include the possibility of atmospheric CO 2 enrichment in our modeling. Third, wildfire is expected to increase regionally across the west due to the expansion of exotic grasses in areas currently in a semi-degraded condition (Creutzberg et al. 2015); however, our projections do not include potential feedbacks related to fire. Relatedly, we assume continuation of current management practices (e.g., grazing intensity) and vegetation treatments. The interactions between grazing and climate observed in the reference period are presumed to continue, for example, proportionally higher grazing intensity in drought years in some locations. Similarly, Creutzberg et al. (2015) assumed the continuation of continued fire suppression and grazing intensity. Fourth, climate change projections have limitations themselves, containing simplifying assumptions and error (Collins et al. 2013). For example, climate models differ in their projections of spatial patterns of precipitation change (Lucht et al. 2006), and there is a high degree of uncertainty in estimates of summer drought (Creutzberg et al. 2015).
Modeling each component separately caries the assumption of independence among components; however, this is not the case. For example, environmental effects often result in nonlinear responses to biological data, frequently resulting from interaction with other species (Clark et al. 2020). It is often difficult to incorporate the effects of these interactions into models. We analyzed the cross-correlations (i.e., bare ground vs. herbaceous cover, shrub vs. sagebrush) among the component suite in the BIT reference period and all six sets of projections. We found similar cross-correlations in all sets of data, with correlations found in the scenarios deviating an average of 8-12% from the BIT. This confirms that the GAM predictions reflect the multicollinearity among components equally well to the training data. Our training data (i.e., the BIT) were modelled in such a way that the summation across components was constrained; the pixel sum v www.esajournals.org needed to be >90 and <110. The projections reflect the constraints of the training data and do not show unrealistically high or low summations or component interrelationships.
Finally, one critical assumption is that vegetation will respond to future climate conditions similarly to observed in the reference period and that no change in soils or topography will occur. While little truly novel climate exists in the future (relative to the spatiotemporal variation in the reference period, Fig. 2), it is possible that prolonged (i.e., decades) departure from historical climate conditions would result in differing vegetation response to that of more short-term variation about the long-term mean observed in the reference period. However, recent climate trajectories have generally tracked those of RCP 8.5 (Sanford et al. 2014), suggesting this signal may be captured in our reference period.

Data analysis
We evaluated the relative importance of each independent variable in model development, assessed with chi-squared values. Next, we plotted mean GAM component cover prediction by given values of independent variables when all other variables are held at their means. For each evaluated time period, we calculated the mean and standard deviation of climate variables and component cover. We plotted the trends observed in the 1985-2018 reference period and the 2020-2080 projection period for RCP 4.5 and 8.5 scenarios. Next, we calculated the climate variable and component cover difference from the 1985-2018 reference period. Third, we evaluated the differences in component cover histograms between the reference period and each GAM projection. With this approach, we can evaluate the shifts in component cover from a density perspective. Finally, we calculated the mean component cover by U.S. Environmental Protection Agency (EPA) level III ecoregions (Omernik and Griffith 2014).

Climate trends
Climate conditions were profoundly altered across much of the study area in the future scenarios (Figs. 2, 3, Table 2). While mean GSPRCP and NGSPRCP were nearly flat over the 1985-2080 period (Fig. 3), increasing by 5.1 and 21 mm, respectively, some regions had dramatic changes of over 100 mm relative to the reference period (Fig. 2). Increased GSPRCP by the 2080s was observed in mountainous areas of Wyoming, southern Montana, the Northern Great Basin, and much of Nevada (Appendix S1: Table S1). Increased NGSPRCP was more widespread, especially in higher elevation mountain ranges across the study area. Decreased GSPRCP and NGSPRCP were common in New Mexico, Colorado, and Utah. Temperature variables increased in all locations, in both RCP scenarios. Mean GSTMAX and NGSTMIN in the RCP 8.5, for example, were projected to increase 6.1°and 4.4°C by the 2080s, respectively (Table 2). GSTMAX increase was relatively homogenous, while NGSTMAX tended to increase more in the south than north. GSTMIN and NGSTMIN increases were strongest in the Northern Great Basin and Wyoming (Fig. 2). The transition between the data observed in the reference period and the projections of the future is smooth across all components (i.e., the trend line of the reference period is nearly parallel with those of the RCP scenarios), with the exception of GSTMAX where the trend line for the reference period is relatively flat relative to the sharp increase projected in the future (Fig. 3). The dramatic amount of interannual variation and large spatial standard deviation in the reference period is also notable.

Component cover trends
The climatic changes were associated with projected changes in component cover (Figs. 4-6, Table 3, Appendix S1: Table S2). Shrub cover was projected to increase from 2020 to 2080s, with similar results between RCP scenarios. The decreases observed to sagebrush cover in the reference period were projected to become stronger through the 2080s, a 2.4% decrease (or~35% proportional reduction) was expected with RCP 8.5 by the 2080s. Decreases to sagebrush cover were projected across much of its range, but some increases were projected in Wyoming, the Northern Great Basin, and eastern Montana (Fig. 5, Appendix S1: Table S2). Ostensibly erroneous increases in sagebrush were projected in southern Nevada, northern Arizona, and southeast Utah due to the GAM flattening the response. Herbaceous cover and annual herbaceous cover both increased in the reference period but were projected to have diverging trends through the v www.esajournals.org 2080s, with the former decreasing substantially, and the latter increasing slightly. The largest increases in annual herbaceous cover were projected to occur in Wyoming, central Oregon, and east-central Idaho as an expansion around core areas originally dominated by invasive plant species during the reference period. Litter cover was largely expected to follow the same directions of herbaceous cover. Bare ground was projected to increase substantially in both RCP scenarios, a net 4.6% increase by the 2080s in RCP 8.5. Bare ground increases were especially great in the driest and lowest elevation portions of the study area, in the Colorado Plateau, Central Great Basin, Snake River Plain, and Great Plains (Appendix S1: Table S2). Decreased bare ground was projected in the Rocky Mountains, Wasatch Mountains, and Blue Mountains. Results were mostly consistent between RCP scenarios for shrub and annual herbaceous projections, but differed in sagebrush, herbaceous, litter, and bare ground, where the RCP scenario projected more dramatic trends. As with the climate variables, the projected component cover trends were nearparallel with the trends observed over the reference period (Fig. 4). Component cover predictions and differences from reference period means in RCP 4.5 and 8.5 scenarios are depicted in more detail for each component in Appendix S1: Tables S1 and S2 and Appendix S1: Fig. S2.
Component cover changes in RCP scenarios were highly interrelated in many cases. For instance, bare ground change between the reference period and RCP 8.5 scenario in the 2080s was strongly related spatially to change in herbaceous (r = À0.67), litter (r = À0.70), and shrub cover (r = À0.57) in the same scenario. Pixels with projected increased bare ground cover were then expected also to have decreased herbaceous, litter, or shrub cover, and vice versa. This demonstrates that while the GAMs were developed independently, they successfully capture the historical range of observed component changes.

Model validation
We used the GAMs to predict component cover using the weather data from 1990, 2000, and 2010. The predictions were compared to the respective years of BIT data. We found similar mean component cover predictions between the GAM and BIT data (Table 1, Fig. 4). However, GAM standard deviation values tended to be smaller than BIT, indicating that GAM predictions were somewhat flatter, possibly due to a lack of imagery input or model generalization. Similarly, GAM predictions have a lower range of component cover predictions relative to BIT (Appendix S1: Fig. S2), though these portions of component histograms excluded by GAMs tended to be quite uncommon in the reference period (Fig. 6, gray shading). Landsat imagery provided the key inputs to the BIT data and was critical to increasing the range and detail of predicted values beyond what is possible with the exclusively nonspectral data used in the current study. The spatial correlations between GAM and BIT predictions (i.e., training) are strong at a broad scale, however, averaging a correlation coefficient (r) of 0.63. Herbaceous and bare ground models were the strongest performers (Table 1), likely due to their robust relationships with GSPRCP and soil variables (Fig. 7.)

Component histogram changes
Projected component changes influence different portions of component histograms, for example, increases in abundance of pixels above x% cover and decreased abundance below y%. Shrub cover is projected to have lower amounts of pixels with less than~8% cover and more at greater cover values, especially around 10-20% cover. Relatively dense sagebrush stands (>10% cover) are projected to become less common, especially in the RCP 8.5 scenario, and by the 2050s and beyond. The large decrease in abundance of 0-2% sagebrush cover is a result of the GAMs flattening the range of values. Herbaceous cover is projected to increase at the tails of the (GS) and nongrowing season (NGS) in the 2080s in RCP 8.5 scenario relative to the 1985-2018 reference climate conditions. Purple indicates novel climate conditions in the 2080s in RCP 8.5 scenario by variable, defined as conditions that were never observed in the reference period throughout the sagebrush biome. Note that color ramps are condensed for clearer symbolization and the true range of values is greater. (Fig. 2. Continued) v www.esajournals.org Fig. 3. Mean values of climate variables for reference period (1985-2018, blue line) and projected climate under RCP 4.5 (red line) and RCP 8.5 (black line). Error bars represent AE 1 SD through space in each time period. Coefficients and R 2 of linear models used to fit the yearly means of the reference period, and coefficients for the RCP 4.5 and RCP 8.5 in blue, black, and red texts, respectively. reference period distribution, with increased abundance of cover above~25% and below 5%, having the effect of broadening and flattening the histogram. Annual herbaceous cover is projected to be increasingly dominated by pixels having greater than 5% cover. Bare ground has complex projected change in cover abundance; increase between~20 and 50, and above~65%; and deceased abundance elsewhere. A notable feature of the RCP 8.5 scenario in the 2080s is the large increase in abundance of pixels with over 80% bare ground cover.

Component cover relationships to independent variables
To evaluate the relationships between component cover predictions from the GAMs and independent variables, we evaluated the mean component cover prediction at given values of each independent variable (with all other variables held constant; Fig. 7). For example, the mean bare ground cover when clay 0-30 cm is 0%, 20%, and so on, with all other variables held at their means. These relationships provide insight into the functional dynamics of component cover in relation to biophysical variables, the core of the GAMs. The range of values on the y-axis of each plot is related to the importance of that independent variable, with a higher range indicating greater importance.
Results from this analysis indicate many significant (P < 0.05), often nonlinear, relationships exist between component cover and biophysical/climate variables (Fig. 7). Full interpretation of these relationships is beyond the scope of this paper, but we highlight a few below. For example, deeper pools (31-300 cm) of OM are negatively related to bare ground cover, but positively related to shrub and especially herbaceous components. Climate relationships were generally as expected and indicate very complex interrelations. While GSPRCP was positively related to all vegetation components, vegetation cover tended to decrease above~1000 mm of NGSPRCP. This is related to these conditions of high NGSPRCP occurring as snow at high elevations (also indicated in the DEM plot of Fig. 7), resulting in short growing seasons and little vegetative cover. Another interesting pattern was the inverse parabolic and (parabolic) relationships with GSTMAX in vegetation components (and bare ground) where peak vegetation cover occurs between~15°and 30°C. For bare ground and (vegetation components), the highest and (lowest) cover values are associated with dry GSPRCP, wet NGSPRCP, low/high GSTMAX, low NGSTMAX, low GSTMIN, and high NGSTMIN. Topographic variables aspect and position index (POSINDEX) were weakly related to component cover, but DEM and slope were both positively related to bare ground cover.

Importance of independent variables to model development
Based on the chi-squared statistics, elevation was the most important overall variable, followed by NGSPRCP, GSPRCP, and OM 0-30 cm (Fig. 8). In general, the 0-30 cm depth soil variables were more influential in model development than the deeper, 31-200 cm variables. Only the climate variables were altered when applying GAMs to future scenarios, but the DEM and soil variables were critical in serving as spatial modifiers because their native resolution was 30 m. 6.8 (3.0) 7.9 (2.9) 8.9 (2.9) 9.6 (2.9) 8.0 (2.9) 9.8 (2.9) 11.9 (2.9) v www.esajournals.org Fig. 4. Mean component cover for the reference period (1985-2018, blue line) and projected cover by generalized additive models (GAMs) under RCP 4.5 (red line) and RCP 8.5 (black line) scenarios. Error bars represent AE 1 SD through space in each time period. The 1990 GAM component cover means (red triangles) and standard deviation are plotted for comparison to 1990 BIT data. Coefficients and R 2 of linear models used to fit the yearly means of the reference period, and coefficients for the RCP 4.5 and RCP 8.5 in blue, black, and red texts, respectively. Specifically, the site-specific conditions of DEM and soil variables could either dampen or enhance the effect of changing climate on a fine scale.

Model development
Drylands constitute ideal environments to test the impacts of climate change in vegetation as they are strongly governed by the spatiotemporal variably in water availability ). Our GAMs were empirically trained on this variability across a vast range of biophysical conditions throughout the sagebrush biome. The responses of component cover to climate variable changes in the spatial and temporal dimensions, and to DEM and soil variables over the spatial dimension are used to develop the GAM for each component. Critical to successfully developing the GAMs is the length of the Landsat/BIT record, which includes multiple wet/drought cycles (Tollerud et al. 2020), resulting in a robust model with minimal extrapolation to future weather. Also critical is the ability of the BIT data to capture gradual change related to interannual weather variation, which was observed iñ 77% of pixels (Rigge et al. 2021). The GAMs effectively captured this variation into empirical models (Table 1). Indeed, cross-validation of GAM results for 1990 relative to BIT predictions yielded robust relationships, though a tendency of the GAMs to flatten the range of component cover must be noted (Table 1). Our model parametrization does not consider the influence of known disturbances in the reference period, and therefore, these effects are also excluded in the future. Thus, our projected component cover likely inflates vegetation cover overall, especially shrub and sagebrush cover, relative to expectations of increased burned area in the future (Creutzberg et al. 2015). Future work to integrate fire influence would be beneficial because increased biomass (that may more quickly dry) is a recipe for increased fire (Creutzberg et al. 2015, Westerling 2016, Tietjen et al. 2017).

RCP scenario comparison
We found materially similar component cover trends between RCP 4.5 and 8.5 scenarios (Fig. 4, Table 3, Appendix S1: Fig. S2), despite large departures in projected climate conditions (Fig. 3, Table 2). In other words, the differences in projected climate by scenario (Fig. 3) were proportionally greater than those of projected component cover (Fig. 4). Similar results were Table 3. Component cover means (and standard deviation) from rangeland areas of the sagebrush biome for the reference period of 1985-2018, and representative concentration pathway (RCP) 4.5 and 8.5 scenarios for the 2020s, 2050s, and 2080s.  reported by Tietjen et al. (2017) who found RCP 4.5 and 8.5 scenarios gave similar results globally, and by Soulard and Rigge (2020) in the northern Great Basin. Our data do show more dramatic changes in the RCP 8.5 scenario, however, with greater increases in shrub and bare ground, and more abrupt decreases in herbaceous cover and sagebrush cover (Fig. 4, Table 3). Moreover, the RCP 8.5 scenario results in more frequent and extreme droughts (Ford et al. 2019), which may not be directly manifest in our results because projected climate is based on means. Historical emissions have closely tracked RCP 8.5 and are expected to match this scenario through the midcentury (Schwalm et al. 2020).

Overall patterns
Widespread trends in climate were ongoing in the reference period (Fig. 2), and across longer timespans (Gonzalez et al. 2010). Past changes are a risk factor for future vulnerability because there is often a time lag between shifts in climate and community structure (Gonzalez et al. 2010). Myriad positive and negative feedbacks exist between altered vegetation community composition and climate (Harte et al. 2015, Tietjen et al. 2017; for example, grass-dominated areas in sagebrush ecosystems intercepted less precipitation, had less transpiration, and increased soil evaporation relative to shrub-dominated areas (Bradford et al. 2014). Even in regions with expected increased precipitation, vegetation response to climate change may worsen ecological drought. However, rangeland vegetation has been subject to intense selection pressure for traits giving resistance and reliance to drought/ aridity and grazing (Quiroga et al. 2010). In summary, vegetation response to climate change is nuanced, with changes occurring at patch scale impacting landscape-scale patterns, and vice versa. Fractional component data are conducive to resolving most of these changes (e.g., densification of a shrubland patch), parameterize models to climate data, and to make meaningful management decisions.
From a broad perspective, our results point to more shrub, bare ground, and annual herbaceous cover and less sagebrush, herbaceous, and litter cover. The projected shift of increasing shrub cover replacement of semiarid grasslands is likely to increase ecological drought (Wilson et al. 2018) because shrubs tend to use more water than herbs. The simultaneous increase in shrub cover and decrease in sagebrush cover indicate an increasing dominance by non-  . Similarly, the projected slight increase in annual herbaceous cover and decrease in herbaceous cover overall indicate an increasing annual to perennial grass ratio, driven by increased climate suitability for invasive annuals. Those regions currently outside the current core areas impacted by invasive grasses are expected to become more climatically suitable, namely the Northwestern Great Plains, Northwestern Glaciated Plains, and Arizona/New Mexico Mountains (Appendix S1: Table S2).

Component specific results
While it is simple to assume that plant communities will shift to higher latitudes and elevation in response to climate change, the true response may be more complex, and vary based on landscape-level hydrological and soil gradient (Debinski et al. 2010, Munson 2013, Munson et al. 2016b. For example, it appears the downside for NPP in drought is greater than the upside in wet periods; in Idaho sagebrush steppe, irrigated treatment plots with 50% more growing season moisture than control plots, NPP increased by 15%, while NPP was reduced by 25% in drought treatment plots with 50% less moisture than control (Tredennick et al. 2018). Furthermore, it is important to note that variation in component cover may not always relate to change in community composition but could instead be related to variation in cover within a species among years, in addition to actual compositional change . Below, we discuss the changes projected for each component group in the context of expectations reported in the literature.
Shrub and sagebrush.-Woody vegetation has been increasing in many temperate drylands globally (Tietjen et al. 2017), in historical grasslands of the southwestern United States (Van Auken 2000), and our future projections. This increase has been linked with chronic excessive grazing (Van Auken 2000). However, shrub cover has also been found to generally increase under drought stress and warming (Debinski et al. 2010, Harte et al. 2015, Gremer et al. 2018).
At the same time, precipitation amount, timing, and type are critical to shrub cover abundance, with winter-dominated precipitation regimes favoring shrubs over grasses (Cook andIrwin 1992, Paruelo andLauenroth 1996). The climate data show that the NGSPRCP proportion of total precipitation increases somewhat through the 2080s (Fig. 3), possibly tempering the influence of increased temperature on sagebrush cover declines.
Sagebrush may be successful in occupying a broad range of semiarid conditions because it has evolved strategies to buffer against precipitation variability and facilitate competitive advantages over herbs (West and Young 2000, Debinski et al., 2010, Tredennick et al. 2018. Indeed, sagebrush was found to be resilient to small changes (+0.2°C) in temperature (Renwick et al. 2017). However, in some regions, these adaptations may prove insufficient to allow the continued propagation and seedling survival of sagebrush (Shriver et al. 2018) in the altered climate conditions of the future. Many locations, especially in the southern edge of the sagebrush extent, will shift from winter precipitation occurring as snow to predominately rain (Klos et al. 2014) to the detriment of sagebrush seedlings (Palmquist et al. 2016b).
Our data demonstrate a clear divergence in expected outcomes for shrub and sagebrush (Figs. 4,6), driven by their diverging responses to climate (Fig. 7). While functionally similar, and having similar responses to GSTMAX, GSPRCP, and NGSPRCP, shrub and sagebrush have unique responses to NGSTMIN, NGSTMAX, and GSTMIN (Fig. 7). NGSTMAX in particular is positively related to shrub cover, with largely the opposite impact to sagebrush cover. However, in cool GSTMAX sites (<18°C), shrub and sagebrush cover responded positively to warming, but in warm sites (>25°C) they responded negatively. Additionally, in cool NGTSMIN sites (<0°C), warming is beneficial to sagebrush, but not shrub, growth. This finding is similar to Kleinhesselink and Adler (2018) who reported that warmer than mean conditions in cool locations increased sagebrush growth and Tredennick et al. (2016) who found warming to benefit sagebrush populations in a relatively cool Wyoming study area. The net result is a 2.4% decrease (or 35% proportional reduction) to sagebrush and a 0.9% increase (or 5.8% proportional increase) to shrub cover RCP 8.5 by the 2080s, indicating an increasing proportion of non-sagebrush shrubs (creosote, rabbitbrush, mesquite, bitterbrush, etc.). This replacement of sagebrush with the aforementioned non-sagebrush shrubs is commensurate with a biome shift from cold desert to warm desert (Sonoran and Mojave Desert), suggesting a northward movement of the ecotone by the 2080s.
Sagebrush responses to climate change vary regionally based on antecedent conditions and projected change. In central Oregon, Creutzberg et al. (2015) found an increase in the extent of moist sagebrush steppe and decrease in dry shrub steppe under all scenarios through 2100. While in southwest Wyoming, shrub and sagebrush cover was projected to decrease by 2050 (Homer et al. 2015). All four models used by Renwick et al. (2017) found that >75% of sites would have increased sagebrush performance in the future. Other investigators reported a more pessimistic outlook; Still and Richardson (2015), for example, reported a 39% loss in suitable climate for sagebrush by 2050. Our sagebrush projections are within the bounds of previous work, with projected increased cover in only 37% of pixels. Our data indicate increased sagebrush cover over the historically cooler and/wetter portions of its range (Appendix S1: Table S2), following the expectations of previous modeling studies (Renwick et al. 2017, Flerchinger et al. 2019. Herbaceous and annual herbaceous. -Herbaceous and annual herbaceous components responded similarly to climate, but the latter was benefited by increased NGSTMIN, and the former sharply decreased at high NGSTMAX (Fig. 7). In arid and semiarid landscapes of the western United States, the distribution of exotic annual grasses has been shown to covary with precipitation patterns at plot and landscape scales measurably more so than native vegetation when temperatures are favorable (Chambers et al. 2007, Pilliod et al. 2017). However, little evidence of this pattern exists in our data, with similar responses to NGSPRCP and GSPRCP (Fig. 7).
While climate projections for the region call for a warmer and slightly wetter climate, cheatgrass cover was predicted to remain stable in >80% of the Northern Great Basin region (Boyte et al. 2016). However, our RCP 8.5 scenario data show increased annual herbaceous cover in 62, 57, and 58% of the study area in the 2020s, 2050s, and 2080s, respectively. Our future projections of annual herbaceous cover are based on climate relationships alone, and many areas are still susceptible to initial invasion and densification and to increased cover following expected increases in burned area. Thus, the projections may underestimate future annual herbaceous cover. However, some of the warmest margins of the sagebrush biome may become more resistant to cheatgrass invasions in the future (Bradley 2009, Creutzberg et al. 2015. Instead, with increased temperature and moisture, these areas may become more suitable to another invasive grass, red brome (Bromus rubens) (Salo 2005).
Bare ground and overall productivity.-Bare ground is the inverse of all the vegetative cover components. In this manner bare ground serves as a convenient shorthand inverse of NPP, the projected increase in bare ground strongly suggests lowered NPP in the future. This inverse relationship can be readily observed in the component change projections (Fig. 5) and in the relationships with biophysical and climate variables (Fig. 7). The projected small net increase in NGSPRCP would be generally expected to decrease bare ground cover, based on the relationships in Fig. 7. The increased temperature variables have diverging influences on bare ground, with negative associations at low GSTMAX and with NGSTMAX (i.e., at mesic sites), but generally positive associations. In general, these patterns follow expectations as more mesic sites are less strongly limited by moisture availability, and instead limited by interspecies competition and/or temperature (Munson 2013, Chen et al. 2019. Indeed, the precipitation to NPP relationship is stronger in the western, drier Great Plains compared with the wetter, eastern Great Plains (Chen et al. 2019).

CONCLUSIONS
Our GAMs capture the empirical spatiotemporal variability of component cover to weather conditions, given the biophysical context of soils and topographical data. The vast range of the sagebrush steppe, variable weather conditions, and myriad biophysical conditions creates a v www.esajournals.org robust GAM for each component, with trivial extrapolation to future climate conditions, and strong cross-validation results relative to BIT data. Our 30-m projections of component cover represent a major improvement upon previous modeling efforts in terms of the robustness and comprehensiveness of training data and the spatial resolution and geographic extent of model outputs. GAMs proved to be more effective than RT in our study area and offer the benefit of a more readily interpretable process. Though the RCP scenarios considerably varied in the projections of future conditions, in general they yield similar results. Both point to plant communities increasingly occupied by non-sagebrush shrubs, annual grasses, and bare ground. Our modeled decrease in sagebrush cover is within the range of findings from previous work, with increased cover by the 2080s modeled in 37% and 44% of pixels in the RCP 8.5 and 4.5 scenarios, respectively. However, the continued widespread conversion of sagebrush steppe to annual grass through fire in landscapes with low ecosystem resistance and resilience (Coates et al. 2016) will likely result in further declines than when considering climate alone.
We confirm the dominant role of temperature changes to the spatiotemporal variability of component cover in the reference period, and thus in the projections. Moreover, the magnitude of temperature changes tends to be proportionally greater in both RCP scenarios than that of precipitation. We demonstrate that in cool sites some degree of warming GSTMAX or NGSTMIN could be beneficial to sagebrush and shrub growth. However, warming NGSTMAX was beneficial to shrub but not sagebrush growth.
Our data are subject to several limitations; they do not consider the impact of disturbance, the impact of interannual weather variability in the future, or changes to management practices. Moreover, our models are designed in such a way to avoid spatial or temporal autocorrelation; due to the serial dependence of the component and climate data, capturing these effects could improve predictive performance (e.g., Vaclavik et al. 2011), and these effects could be considered in future work. However, our results follow expected spatial patterns and provide information useful to managers in preparation for future vegetation composition and cover through the prioritization of conservation and restoration and shed light on species range shifts.