Predicting Exotic Annual Grass Abundance in Rangelands of the Western United States Using Various Precipitation Scenarios

ABSTRACT Expansion of exotic annual grass (EAG), such as cheatgrass (Bromus tectorum L.) and medusahead (Taeniatherum caput-medusae [L.] Nevski), could cause irreversible changes to arid and semiarid rangeland ecosystems in the western United States. The distribution and abundance of EAG species are highly affected by weather variables such as temperature and precipitation. The study's goal is to understand how different precipitation scenarios affect EAG abundance estimates and dynamics, and we develop a machine learning modeling approach to predict how changes in annual and immediate past precipitation patterns could affect the abundance of EAG. The machine learning predictive model used seed source from previous years, weather variables, and soil profiles to drive its predictions. We achieved excellent training accuracy (r = 0.95 and median absolute error [MdAE] = 2.36% cover) and strong test accuracy (r = 0.79 and MdAE = 4.54% cover). We developed five versions of EAG abundance maps for 2022 with different precipitation scenarios: 9 yr of average precipitation, half of the average, three-fourths of the average, one and one-half times the average, and two times the average. The approach presented can be replicated to new study domains and easily modified for use with other precipitation scenarios. Developing multiple versions of a year's EAG spatially explicit abundance dataset predictions from multiple weather-based scenarios can provide important information to land managers as they prepare for variable EAG dynamics each year. Informed annual predictions based on weather scenario–driven models have the potential to improve fire preparation decisions.


Introduction
Invasion of exotic plants, such as cheatgrass ( Bromus tectorum L.), red brome ( Bromus rubens L.), and medusahead ( Taeniatherum caput-medusae [L.] Nevski), has led to the degradation of ecological services and reduced biodiversity in arid and semiarid rangelands of the western United States ( Bradley et al. 2017 ;Pilliod et al. 2017 ). Exotic annual grass (EAG) species are dominant invaders in these ecosystems, affecting grazing land and wildlife habitat by altering fire regimes ( D'Antonio and Vitousek 1992 ;Balch et al. 2013 ), which, consequently, reduces the abundance and diversity of native plant species ( Davies et al. 2011 ;Pyšek et al. 2012 ;Tarbox et al. 2022 ). Residual plant materials of dead EAGs provide abundant fine fuels and contribute to increased wildfire frequency and intensity ( Bradley et al. 2017 ), which exacerbates habitat loss. The native habitat loss presents substantial challenges for land managers who work to maintain and restore native systems in these dry environments because habitat restoration in these ecosystems is costly, time consuming, and many times unsuccessful due to the variability of precipitation patterns ( D'Antonio and Vitousek 1992 ; Davies et al. 2011 ;Chambers et al. 2013 ;Bradley et al. 2017 ;Germino et al. 2018 ;Smith et al. 2021 ). Winter annual EAGs, with shallow rooted systems, often start their growing seasons in fall if they receive adequate moisture. If EAGs do not receive adequate moisture during fall, they will germinate and emerge in winter and/or early spring when soil moisture interacts with warmer temperatures and environmental conditions favor growth. The early start of EAGs' growing season gives them an advantage over many native plant https ( Chambers et al. 2013 ;Pilliod et al. 2017 ).
EAG distribution and abundance are strongly affected by topography, such as elevation and hillside position; competition; and precipitation regimes ( Bradley and Mustard 2005 ;Bykova and Sage 2012 ;Chambers et al. 2013 ). Several studies have shown that a changing climate may alter future weather patterns and fire regimes, which could more rapidly increase the spread of EAGs at higher elevations ( Bansal et al. 2014 ;Compagnoni and Adler 2014 ;Smith et al. 2021 ). Rapid climate change can make it difficult to predict what future climate regimes will bring in terms of weather, particularly with regard to precipitation, and how these changed regimes will affect the distribution and abundance of EAGs ( Bykova and Sage 2012 ;Bradley et al. 2016 ;Pilliod et al. 2017 ). An approach to EAG prediction and mapping that considers a range of precipitation scenarios can provide a more holistic view to the management of EAGs and related policy. Although a few laboratory-based studies have been conducted to see what effect different moisture scenarios have on EAG species ( Bansal et al. 2014 ;Bishop et al. 2020 ;Kainrath et al. 2021 ), to our knowledge no study has predicted EAG abundances using multiple precipitation scenarios on a large geographic scale. However, several studies estimated EAG abundance at various temporal and spatial scales Bateman et al. 2020 ;Pastick et al. 2020 ;Pastick et al. 2021 ;Larson and Tuor 2021 ;Weisberg et al. 2021 ;Dahal et al. 2022 ).
In this study, we set out to develop an EAG model prediction method with resulting maps that use multiple precipitation scenarios. The model and 30-m resolution maps cover a large regional scale using an approach that could be replicated to any spatial domain. The goal is to understand how different precipitation scenarios affect EAG abundance estimates and dynamics. The study objectives are to 1) develop a model that can ingest multiple alternative precipitation scenarios and generate associated spatially explicit abundance maps, and 2) analyze the spatially explicit maps to understand how they reflect the dynamics of the different precipitation scenarios, with a focus on effects of elevation. We hypothesize that more precipitation will lead to increased EAG abundance, but that a precipitation threshold may occur in some places where EAG abundance will cease to increase and may even decrease.

Study area
The study area consists of all or part of 17 western states in 33 level III ecoregions ( Fig. 1 ). The ecoregions provide a geographical framework for ecosystem management and environmental understanding in the broadest sense including abiotic and biotic components of terrestrial and aquatic ecosystems ( Omernik and Griffith 2014 ). Elevations range from 86 m to 4421 m ( Gesch et al. 2018 ). The entire study area spans multiple altitudinal zones including plains, foothills, montane, subalpine, and alpine vegetation communities. However, the study only includes pixels classified as grassland/herbaceous or shrub by the National Land Cover Database (NLCD) 2019 ( Jin et al. 2019 ; Dewitz and US Geological Survey 2021 ) that were at or below 2350-m elevation, which are likely arid and semiarid rangelands. The rangeland area is dominated by native shrub species (primarily sagebrush [ Artemisia spp.]) and herbaceous vegetation containing both native and exotic perennial and annual forbs and grasses. Native perennial shrub and bunch grass communities in much of the northern part of the study area green up in early to mid-May and remain photosynthetically active through the hot summer months. The exotic annual grasses can green up earlier in spring and senesce and cure by late spring or early summer ( Brooks et al. 2004 ;Bradley and Mustard 2005 ;Bradley et al. 2016 ;Rigge et al. 2021 ;Dahal et al. 2022 ). Some native species in the southern part of the study, however, might go dormant during the summer months, phenologically coinciding with exotic annual grasses ( Benedict et al. 2023 ;Salo 2004 ). The United States has experienced changes in precipitation and temperature patterns since the beginning of the 20th century. Average annual precipitation in the northern part of the study area has experienced upward trends while downward trends have affected the southern part. These trends are expected to continue ( Hayhoe et al. 2018 ). The 9-yr average (2012-2021) annual precipitation in the ecoregions ranged from 150 mm to 1125 mm. The summer average temperature of the ecoregions ranged from 23 °C to 41 °C, whereas winter average temperature ranged from −13 °C to 4 °C (PRISM Climate Group: https://www.prism.oregonstate.edu ).

Input data and preprocessing
We acquired data from 23,287 field survey plots developed by the Bureau of Land Management (BLM) Assessment, Inventory, and Monitoring (AIM) program (i.e., the Landscape Monitoring Framework [LMF] and Terrestrial Aim Database [TerrADat]) from 2016 to 2021 ( Fig. 1 ) to be used as a regression-tree model's response variable. The AIM program uses standardized monitoring protocol to collect field data with a line-point intercept method, which focuses on measuring core terrestrial indicators including plant species cover and composition for plots chosen using a probabilistic sampling design based on AIM core indicators, such as management concern, soil stability, vegetation composition and height, and bare ground, within federally managed lands. To collect the indicator variables, measurements of 101 pin drops in two 150-ft intersecting transects (LMF) and typically 150 pin drops in three transects (TerrADat) are taken and aggregated for each plot ( Herrick et al. 2005 ;Kachergis and Burnett 2021 ). We compiled 16 individual EAG species from the AIM database and aggregated them to create our model training dataset (for more detail see Dahal et al. 2022 ).
Weather variables at a monthly timestep-precipitation (ppt), minimum temperature, and maximum temperature-were acquired for 2012-2021 from the PRISM Climate Group ( https://www.prism. oregonstate.edu ). Raster layers depicting each monthly weather variable at 800-m resolution were reprojected to NAD83 Albers and resampled to 375-m resolution using the cubic convolution method. The layers were again resampled to 30-m resolution using the same method to align with other datasets used in the model and subsequent analysis. We used a multistep resampling approach to minimize hard lines in output products that can be evident when inputs of different native resolutions are used in model development ( Cracknell 2010 ;Gomez-Chova et al. 2011 ).
From the monthly weather variables, we calculated several derivative layers that were input into the regression-tree model and used as predictor variables in mapping the final outputs. Derivatives included the annual water year's (previous Octobercurrent September) monthly ppt, maximum temperature, and minimum temperature layers. We also calculated the long-term (2012-2021) average annual and average monthly layers from the monthly weather variables. Following Pilliod et al. (2017) , we calculated 2 previous yrs' average (2yrPAve) and 3 previous yrs' maximum (3yrPMax) ppt layers. To map EAG abundances with variable ppt scenarios, we calculated five ppt layers using the 9-yr average (normal, or 1.0) and factors of that 9-yr average including 0.5, 0.75, 1.5, and 2.0, where, for example, a factor of 0.5 (ppt_0.5) inputs values half of the normal ppt and a factor of 2.0 (ppt_2.0) inputs values double the normal ppt.
Along with weather data, we input several other variables to drive the predictive model. A previous year's EAG abundance data were used as proxy for a dynamic seed source layer ( Dahal et al. 2022 ). Topographic variables such as elevation ( Gesch et al. 2018 ) and potential annual incident direct radiation (PADR) ( McCune and Keon 2002 ) were included because EAG performance and distribution are strongly affected by location ( Chambers et al. 2013 ). Bansal et al. (2014) concluded variations in soil properties made noticeable differences on some of the EAG growth, so we input data from five edaphic variables into the model-clay, sand, silt, soil organic matter, and available water capacity-within a 0 to 30-cm depth zone ( Chaney et al. 2016 ).

Modeling and mapping
The model database was created by spatially aligning predictor variables for each AIM field survey plot, coincident with the year of the survey. For example, if an AIM field plot was surveyed in 2020, then the predictor variables(i.e., weather variables, EAG abundance, 2yrPave, and 3yrPmax) would be applicable to 2020. Some variables like elevation and edaphic variables, were static, so they did not change on the basis of the year of the AIM field plot survey data. These static variables were applicable to all years in the study period (i.e., 2016-2021).
The entire database was randomly divided into 80% training data and 20% independent test data for model building and validation purposes, respectively, to make sure the test dataset is representative of the entire population but model robustness and improvement are still achieved ( Wylie et al. 2019 ). A gradientboosting regression-tree model was developed using python scikitlearn and lightGBM software libraries ( Pedregosa et al. 2011 ;Ke et al. 2017 ). We used an early-stop approach and the grid search hyperparameter optimization method to prevent overfitting and identify optimal parameters. To gain different insights from the data, we optimized learning rate, maximum depth, subsample (record sample), and feature-fraction (feature sample) parameters. Optimizing subsample/feature-fraction parameters forced the algorithm to learn from different sets of records/features in each iteration, limiting the overusage of records/features. The "linear tree" parameter was activated while developing the model so that the regression tree could extrapolate beyond the learned values. The calibrated model predictions were compared with the independent test data to validate model performance, and model accuracy results are reported.
To develop the weather-based scenario maps, the calibrated model algorithms were applied to spatial versions of static variables (edaphic, elevation, and PADR) and dynamic variables appropriate for 2022 (monthly and annual weather variables, 2021 EAG abundance, 2yrPAve from 2020 to 2021, and 3yrPMax from 2019 to 2021) using the lightGBM, gdal, and numpy libraries in a Python environment ( https://gdal.org/ and https://numpy.org/doc/ stable/index.html ). For the ppt_1.0 (normal) scenario EAG map, we applied the model algorithms along with the normal ppt data layers and other static and dynamic variables. To develop the other ppt scenarios (e.g., ppt_0.5), we substituted the ppt data associated with a specific scenario for the normal (ppt_1.0) scenario in the model. While mapping, we applied a mask to areas that were not classified as grassland/herbaceous or shrub by the 2019 NLCD and areas above 2350-m elevation to target only likely rangeland ecosystems and areas that are more likely ecologically suitable for EAG species.

EAG in relation to elevation and precipitation scenarios
Several studies have indicated elevation can be an important driver of EAG distribution and abundance ( Chambers et al. 2013 ;Compagnoni and Adler 2014 ;Smith et al. 2021 ). To analyze EAG abundance along various elevation gradients within the context of ppt scenarios, we categorized elevation into five classes: ≤ 10 0 0 m, > 10 0 0 −≤ 150 0 m, > 150 0 −≤ 1750 m, > 1750 −≤ 20 0 0 −m, and > 20 0 0 −≤ 2350 m. We then extracted predicted pixel values of EAG abundance from each ppt scenario to the elevation classes.

Statistical analysis
We evaluated the model performance by calculating Pearson's correlation coefficient ( r ) and median absolute error (MdAE) for training and test samples. To present the results, we generated scatterplots of observed and predicted values of training and test samples (see section 3.1). To understand the importance of each independent variable used in the modeling, we computed partial dependency of the variables on the dependent variable and generated plots for highly used independent variables (see section 3.1). Maps were evaluated in various ways. First, we calculated mean, median, maximum, and standard deviations for each scenario map (see section 3.2). Second, we categorized the ppt_1.0 maps by deciles and evaluated each scenario map by calculating the mean percent values to each decile (see section 3.2). Third, we evaluated the relation of EAG abundance to elevation cohorts and ppt scenarios by calculating mean, median, standard deviation statistics (see section 3.3). Our hypothesis was to see whether ppt scenarios play significant roles in EAG abundance; therefore, we conducted statistical analysis using the Mann-Whitney U test, where we compared the baseline (ppt_1.0) distribution of EAG abundance to other scenarios as a function of elevation. The Mann-Whitney U test is a nonparametric test, which is often used as a test of difference in location between distributions. We generated violin plots for each scenario and elevation class to highlight the full distribution of the data (see section 3.3). We used matplotlib, scikit-learn, and scipy Python modules for the statistical analyses and graphics.

Modeling and validation
Correlation analysis of the predicted values with the independent test samples ( n = 4657) showed strong agreement with an r of 0.79 ( P value < 0.001) and MdAE of 4.54% cover. The r and MdAE for predicted values with the training samples ( n = 18,630) are 0.95 ( P value < 0.001) and 2.36% cover, respectively ( Fig. 2 ). In the test scatterplot, the 1:1 line crosses the regression line at a value of about 15% cover. Below that level of predicted value, the model generally overpredicts, especially when the observed values are close to zero. The model generally underpredicts when observed values are above where the two lines intersect. Far fewer points are underpredicted ( n = 1601) because areas invaded with high abundance are less common than areas invaded with low abundance.
Model agnostic interpretation methods, as partial dependence plots, reveal the feature importance, or how much an independent Figure 3. Partial dependence plots (PDP) from the lightGBM regression-tree model for the six most highly used variables. The plots are arranged in order of decreasing importance from left to right and top to bottom. PreYrEAG is the previous year's exotic annual grass abundance as a proxy for seed source, ppt_pYear_m10 is precipitation of October from the previous year, PADR is potential annual incident direct radiation, and ppt_cYear_m03 is March precipitation for the current year. Silt and sand are soil texture within 0-30 cm from the surface. The hash marks at the base of each PDP represent the deciles of the predictor variable distribution.
variable was used in model development, and unique interactions between the predicted EAG abundance and predictor variables (Fig.  3). Model predictions of EAG abundance varied largely as a function of the previous year's EAG cover (proxy for seedbank and invasion envelope), soil types (e.g., silt and sand), and the previous year's October ppt along with PADR and March ppt of a current year. The relation between the prediction of EAG abundance and previous years of EAG cover was close to sigmoidal with low and high invasion in areas with low and high variability, respectively. Other variables were found to be useful predictors of EAG abundance but, by themselves, were not highly useful for separating each component. These variables had generally similar predictive thresholds around 15% of EAG cover, which is closer to the value where the model showed a higher agreement (see Fig. 2 ).

Scenario mapping
Five ppt-based scenario EAG abundance predictions for 2022 are displayed as maps ( Fig. 4 ), with abundance presented as percent cover ranging from 0 to 100. In general, the most abundantly invaded pixels occurred in the northern tier of the study area and on the perimeter of the Mediterranean California ecoregion (ecoregion 11.1.2 on Fig. 1 ). The least abundantly invaded pixels occurred in the southern tier of the study area, especially on the eastern side. Spatially, as ppt deviates from the normal (long-term average), the prediction of EAG abundance generally follows. Overall, EAG abundance has a positive correlation with ppt: i.e., when ppt decreases, the predicted EAG abundance decreases, and vice versa. This is demonstrated in Table 1 , where the overall mean for the Table 1 Overall (top panel) and nonzero pixels (bottom panel) mean, median, maximum, and standard deviation (STD) of exotic annual grass percent cover by different precipitation scenarios. ppt_1.0 EAG abundance map for the study area was 8.75%. When scenarios for ppt_0.5 and ppt_0.75 were mapped, predictions of the overall mean decreased to 6.86% and 7.64% cover, respectively, relatively substantial decreases. Conversely, when the ppt_1.5 and ppt_2.0 scenarios were mapped, predictions of the overall mean increased to 9.21% and 9.84% cover, respectively, approaching a critical ecological threshold of 10% cover for increased wildfire activity in sagebrush ecosystems ( Pastick et al. 2021 ). Approximately onethird (ranging from 29%-33%) of the study area has zero estimated EAG abundance regardless of the ppt scenario, which reduces the overall mean and median EAG abundance for each scenario. This also indicates that about one-third of the study area is still free from invasion of the 16 EAG species ( Dahal et al. 2022 ) aggregated to the study. If pixels with zero values are removed from the calculation of the statistics, then a more accurate effect of ppt vari- ability on EAG abundance in invaded areas is gleaned. The recalculated statistics with these zero pixels removed from the analysis are shown on the bottom half of Table 1 . Maximum and standard deviation values for the study area fluctuated from 81% to 100% cover and 9.09 to 13.20, respectively, following the same trend as the mean and the median. Lower mean, median, and standard deviation values were associated with lower ppt scenarios, and higher mean, median, and standard deviation values were associated with higher ppt scenarios. The overall median EAG abundance values for all ppt scenarios were similar (range of 3-7) no matter what level of ppt. However, the median values were substantially higher (range of 7-10) when the zero value pixels were removed from the calculation.
Comparison of average EAG abundance of all scenarios within the percent cover deciles of ppt_1.0 (normal ppt) shows EAG abundance decreases substantially more with lower-than-normal ppt scenarios than it increases with higher-than-normal ppt scenarios ( Fig. 5 ). The average ppt_1.0 EAG abundance value within the 90-99% decile is 92.86%, whereas the average values for the ppt_1.5 and ppt_2.0 scenarios are 96.20% and 97.69%, respectively, and the average values are 71.09% and 84.66% for the ppt_0.5 and ppt_0.75 scenarios, respectively. Similarly, the average EAG abundance value for the ppt_1.0 scenario within the 60-69% decile of ppt_1.0 areas is 63.42%, and the values for the same areas of the ppt_1.5 and ppt_2.0 scenarios are 63.72% and 66.05%, respectively, and conversely 46.47% and 56.69% for the ppt_0.5 and ppt_0.75 scenarios, respectively. A total of 2302 pixels were predicted to have > 80% EAG abundance in the ppt_1.0 scenario. The number of pixels with > 80% abundance for the other ppt scenarios equal 1, 76, 31,642, and 194,811 for ppt_0.5, ppt_0.75, ppt_1.5, and ppt_2.0, respectively.

Relation of EAG abundance to elevation and precipitation scenarios
The EAG abundance means increased at all elevation levels across all increasing ppt scenarios except at the > 20 0 0-≤ 2350-m level, where EAG means fluctuated from scenario to scenario ( Table  2 ). For example, mean values at the ≤ 10 0 0-m level ranged from 10.16 to 15.97, steadily increasing as ppt values increased, whereas at the highest elevation level, mean values fluctuated from 3.15, 3.04, 3.40, 3.21, to 3.46 across increasing ppt scenarios. The highest mean values for the lowest four ppt scenarios (ppt_0.5-ppt_1.5) all occurred at the > 1500-≤ 1750-m level. Only in the ppt_2.0 scenario did the highest mean occur outside of the > 1500-≤ 1750-m cohort group, and then it occurred at the < 10 0 0-m level. By a substantial amount, the > 20 0 0-≤ 2350-m elevation cohort group had the lowest means across all ppt scenarios. The EAG medians either remained flat or increased at all elevation levels except, again, at the highest elevation level where the median values fluctuated ( Fig. 6 ). The rank of EAG abundance variability across elevation levels and ppt scenarios varied only between the > 1500-≤ 1750 and the > 1750-≤ 20 0 0 levels, where they each were either the most variable or the second most variable cohort group in all ppt scenarios. Fig. 6 shows the distribution of EAG abundances as a function of elevation and ppt. From the Mann-Whitney U test, we found that the distribution of EAG abundance in the ppt_1.0 (baseline) scenario is significantly ( P value < 0.0 0 01) different than all other scenarios within each elevation class.

Discussion
Building a robust regression-tree model that is balanced to a wide range of conditions is important for developing accurate predictions, but it is challenging because the model needs to capture extreme values and a proportional distribution of the prediction range of the training data and overall population ( Wylie et al. 2019 ). Predictive regression tree models are complex. They are designed to capture higher-order interactions between the input variables for robust modeling and mapping accuracies ( Michaelsen et al. 1994 ;De'ath and Fabricius 20 0 0 ;Ghimire et al. 2012 ). Regression tree modeling, especially with boosted ensemble algorithms, determines, through hyperparameter selection and tree pruning, how to stratify the data into rules and reduce noise to permit the algorithm to use variables when and where they find optimal prediction accuracies ( Wylie et al. 2019 ). We experimented with this strategy to stratify input data by records (subsample) and by columns (colsample) in a lightGBM framework, and, with judicious pruning of regression trees, to find the optimal model- ing parameters. At the same time, we developed the model to reduce the likelihood of overfitting through excessive usage of noisy observations and any individual variable. The resulting modeling accuracies (see Fig. 2 ) and variable usage (see Fig. 3 ) by the optimized model boosted our confidence for developing the EAG abundance maps with different ppt scenarios. The application of a lightGBM modeling framework with linear-tree function enabled the boosted-tree model to extrapolate prediction values beyond what it had learned from training values during the model development. This application produced maps more representative of scenarios ppt_1.5 and ppt_2.0. Without the extrapolation function, the model parameters would restrict the EAG abundance predictions to the training data range and not reflect conditions that the more ample ppt scenarios produce ( Armstrong 2001 ;Gu et al. 2012 ).
Some variables are more important than others when building predictive models. In our EAG abundance prediction model, the previous year's EAG abundance maps, which serve as a proxy for seed source in previously invaded areas, were highly influential variables (see Fig. 3 ). This finding aligned with other studies that noted that once an EAG species gains a foothold in an area, it typically remains present in that area ( Brooks et al. 2004 ;Figure 5. Average exotic annual grass values of all precipitation scenarios in each decile of the normal precipitation scenario. Ppt_0.5 is half of normal, ppt_0.75 is three-fourths of normal, ppt_1.5 is one and one-half times normal, and ppt_2.0 is double the normal (ppt_1.0) precipitation.  Bradley and Mustard 2006 ;Pastick et al. 2021 ). This persistency can be attributed, in part, to the positive feedback cycle between EAGs and fire. Nelson et al. (2014) and Barnard et al. (2019) also noted that soil profiles contribute substantially to the growth and structure of vegetation in the arid and semiarid ecosystems through moisture and nutrient availability. Specifically, topsoil profiles can be important factors in the development and dominance of shallow-rooted EAG species in arid and semiarid ecosystems ( Pastick et al. 2021 ). Bansal et al. (2014) noted edaphic variables played important roles in the growth of three selected EAG species, which aligned with our study's findings where soil profile characteristics (silt and sand) were found to be important variables in predicting EAG abundance.
Another influential variable in our predictive model is the total monthly ppt from a previous year's October, which was the third most highly used variable. Our model shows that increased levels of previous year's October ppt returns higher EAG abundance in the subsequent growing season, which was in line with the finding of Bishop et al. (2020) and Horn et al. (2017) , who noticed that fall ppt significantly increased EAG densities, heights, and biomass by extending the growing season and providing increased germination success. Bentley and Talbot (1951) observed that growth and abundance of winter annuals, like cheatgrass, depend on the ppt distribution during critical growth periods such as fall and spring more than on total annual ppt. March ppt was the sixth most highly used variable in our model. Cheatgrass, which is the most ubiquitous EAG species in our study area ( Pastick et al. 2021 ;Dahal et al. 2022 ), was found to be extremely sensitive to moisture and temperature ( Bradley and Mustard 2006 ;Concilio et al. 2013 ). Cheatgrass's tolerance to upper and lower temperature thresholds guides seed germination, growth, survival, and fecundity. PADR, a proxy for solar radiation, incorporates aspect, slope, and latitude and is another highly influential variable in our model ( McCune and Keon 2002 ;Pastick et al. 2021 ). Solar radiation, along with ppt and elevation, affect soil temperature and water availability, thereby affecting EAG abundance at the landscape level ( Williamson et al. 2019 ). Our predictive model is strongly affected by these factors (solar radiation, edaphic, and ppt), elevating our confidence in map predictions. Compagnoni and Adler (2014) reported increased cheatgrass densities in all their experimental sites that received abnormally high ppt, and the increase was pronounced at higher elevations. We found elevated EAG abundance in some outlier areas in the > 1750-≤ 20 0 0-m elevation class with the ppt_1.5 and ppt_2.0 scenarios (see Fig. 6 ). However, we did not see an overall EAG abundance increase in the > 1750-≤ 20 0 0-m elevation class nor the elevation class above that. When we further analyzed the predicted ppt scenario maps along with the elevation gradients, we found that EAG abundance was higher in low-elevation areas (up to 1750 m) regardless of the ppt scenario. As noted earlier, ppt is the most important variable to EAG abundance predictive modeling because small differences in ppt amounts and distribution can result in large differences in EAG abundance (see Fig. 3 ). We found that lower-than-normal ppt conditions have a larger effect on EAG abundance than do higher-than-normal ppt conditions (see Fig. 5 ). Young et al. (1969) noted that lower-than-normal ppt in fall and spring resulted in greatly reduced germination and elevated mortality of cheatgrass. In an experimental study in Oregon, a 17-fold difference in cheatgrass productivity was recorded between normal and drought years ( Ganskopp and Bedell 1979 ). Bradley and Mustard (2005) noted cheatgrass productivity and abundance could be 10 times higher in a wet year than in a dry year. And although EAG abundance increases to a certain degree when ppt increases, the increase is lower in historically higher abundance areas when ppt is at least one and one-half times or double the average ppt level. Excessive moisture resulting from abnormally high ppt may negatively affect the shallow root systems of EAG species in historically highly invaded areas ( Peek et al. 2005 ;Ryel et al. 2010 ). Conditions such as nitrogen deprivation, seed runoff, and oxygen deprivation could result in less growth and impair the invasion of the shallowrooted grass species ( Tilman and Wedin 1991 ;Garrison and Stier 2010 ).
We noticed historically high EAG abundance areas have higher variation as ppt deviates from normal. A formal validation of the spatial scenario maps remains problematic as no such historical data are available; however, we performed semistructured comparisons of the maps to see whether the model performed as expected considering the ppt scenario (see Table 1 , Figs. 4, 5 ). Our resulting maps agreed with other studies that have documented that EAG abundance follows the amount and intensity of recent past months and years of ppt ( Bradley and Mustard 2005 ;Bykova and Sage 2012 ;Pilliod et al. 2017 ;Bishop et al. 2020 ;Applestein and Germino 2022 ). A visual inspection of the spatial patterns was performed as human eyes can detect spatial patterns that statistical procedures might miss ( Sohl et al. 2012 ). One important piece of scenario modeling is to capture unseen conditions with some sets of defined conditions; however, we also understand that no single approach is best to define the conditions.

Implication
The ecological effect of EAG to the arid and semiarid rangeland ecosystems of the western United States is enormous. This study offers an improved understanding of EAG dynamics with the development of five spatially explicit EAG maps for 2022 using five different ppt scenarios across rangelands of the western United States (see Fig. 4 ). Analysis of the datasets can stimulate spatial understanding useful for land management, such as the ppt threshold where the model estimates that an EAG abundance threshold will be crossed, the 10% wildfire threshold identified by Pastick et al. (2021) , or which increased ppt scenario(s) causes the model to estimate a decrease in EAG abundance. The knowledge gained from this information mitigates the EAG knowledge gap and can help land managers better understand the landscape dynamics of western rangelands and combine it with other information to make more informed land management decisions.