Extracting exotic annual grass phenology and climate relations in western U.S. rangeland ecoregions

This research builds upon the extensive body of work to model exotic annual grass (EAG) characteristics and invasion. EAGs increase wildland fire risk and intensifies wildland fire behavior in western U.S. rangelands. Therefore, understanding characteristics of EAG growth increases understanding of its dynamics and can inform rangeland management decisions. To better understand EAG phenology and spatial distribution, monthly weather (precipitation, minimum and maximum temperature) variables were analyzed for 24 level III ecoregions. This research characterizes EAGs’ phenology identified by a normalized difference vegetation index (NDVI) threshold-based interpolation technique. An EAG phenology metric model was used to estimate a growing season dynamic for the years 2017–2021 for shrub and herbaceous land cover types in the western conterminous United States (66% of the area). The EAG phenology metrics include six growing season metrics such as start of season time, end of season time, and time of maximum NDVI during the growing season. The models’ cross validation results for Pearson’s r ranged from 0.88 to 0.95. Increased understanding of the effects that weather conditions have on EAG growth and spatial distribution can help land managers develop time-sensitive plans to protect entities deemed valuable to society like native habitat, wildlife, recreational areas, and air quality.


Introduction
Monitoring exotic annual grasses (EAGs) using satellite imagery helps researchers gather information at a broad scale to help identify and better understand invaded areas that are at risk of potential irreversible change. EAG species invade various vegetation types including native perennial grasses that compete for the same resources. Some exotic annual grasses such as cheatgrass (Bromus tectorum), medusahead (Taeniatherum caput-medusae), and wiregrass (Ventenata dubia), with shorter root systems, germinate in the winter or early spring taking nutrients from the subsurface (Thill et al. 1984;Weisberg et al. 2021;Wilsey et al. 2017). Native perennial grasses with deeper root systems go dormant in winter and can grow back in spring, after EAGs, when enough moisture is available. EAG emergence dates vary across our study area, influenced by temperature, precipitation timing, ground disturbance, and species (Boyte et al. 2018Bradley and Mustard 2008). Beginning in 2016, the U.S. Geological Survey (USGS) produces annual fractional cover estimates for cheatgrass, medusahead, and 16 EAG species in the western U.S. rangelands (Dahal et al. 2021). As of the 2020 estimates, areas with 15% or more EAG cover occupy over 100 million acres in the western U.S. rangelands (Dahal et al. 2021).
Much EAG research has focused on local environments; however, some researchers have expanded their scope to larger study areas (Bart et al. 2017;Bradley and Mustard 2008;Taylor et al. 2021;Weisberg et al. 2021). These studies have been useful for understanding regional EAG characteristics and phenology. For example, Weisberg et al. (2021) identified medusahead and cheatgrass phenology with high spatial resolution in the Cold Desert ecoregion in northern Nevada using an unoccupied aerial vehicle (UAV). They used the UAV to gather in situ observations during the growing seasons of cheatgrass and medusahead. They noted that the sustained growth during the dates of their imagery occurred between spring thaw and summer dry down, May 15-July 20, 2017. During the selected observation dates, they noted that cheatgrass and medusahead experienced different phenology stages. During cheatgrass active onset, medusahead was small and indistinguishable, but as cheatgrass plants senesced, medusahead began to mature (Weisberg et al. 2021). The observations in the Cold Desert ecoregion can inform our study on estimating sustained growing season for a collection of 16 EAG species. For the purposes of this study, we define sustained growth as normalized difference vegetation index (NDVI) values consistently reflecting increased photosynthetic activity for at least five consecutive weeks.
Different approaches can be used to extract the start of sustained growth of EAG, and one approach that can be relatively easily produced is the thresholdbased approach. The threshold-based approach has been utilized by many studies to extract land surface phenology or EAG phenology (Taylor et al. 2021;Wang et al. 2018;Xie et al. 2022). The approach usually incorporates the amplitude (AMP) from a NDVI time series, which is calculated using the start of the increasing NDVI values and the peak NDVI value within the time series. The technique we used to extract thresholds involves manually interpreting NDVI data profiles that were developed from 2016 to 2020 Harmonized Landsat and Sentinel-2 (HLS) data (Claverie et al. 2018). This moderate spatial resolution (30-m) dataset is suitable for monitoring vegetation and developing training datasets of phenology metrics for an EAG phenology model (Claverie et al. 2018;Dahal et al. 2022;Pastick et al. 2020).
Developing a time series of annual spatial datasets of EAG phenology allows researchers and land managers to track the timing and progression of EAG growing-season dynamics year-to-year. Changing phenology patterns can indicate changing environmental conditions, such as progressively warming temperatures that can cause earlier starts to the growing season (Howell et al. 2020) and extend fire seasons. In the western United States, EAGs have invaded shrub ecosystems for more than a century and, in that time, a positive feedback cycle between fire and EAGs has developed. EAGs increase wildfire frequency and intensity with their compact and combustible tinder-like characteristics (Bradley et al. 2018). According to the Congressional Research Service (Hoover and Hanson 2022), in 2020 and 2021, the United States experienced about 59,000 wildfires each year that burned more than 10.1 million acres and 7.1 million acres, respectively. As fires burn through landscapes, EAG species like cheatgrass take advantage of the resulting bare ground and "cheat" their way into their new home by displacing native perennials and taking up nutrients and available soil moisture following fires Bradley 2013;D'Antonio and Vitousek 1992;Davies et al. 2021;Fusco et al. 2019;Horn and St. Clair 2016).
With the number of fires increasing in the western United States (2019 had over 21,000 fires and that figure increased to over 23,000 fires in 2021 Hanson 2020, 2022), EAGs are spreading and creating more opportunities for wildfires to occur Bradley et al. 2018;D'Antonio and Vitousek 1992;Dahal et al. 2022;Davies et al. 2021;Homer et al. 2012). EAG phenology analysis can help land managers identify when senesced and dead EAGs could be at risk of igniting wildfires.
The timing and magnitude of EAG sustained growth varies between and within ecoregions; and understanding these dynamics can help improve land management and help identify rangeland readiness for livestock grazing (Brownsey et al. 2017;Davies et al. 2022Davies et al. , 2021Thill et al. 1984). EAGs may be used as a grazing resource, but it is not a sustained resource for livestock during an entire growing season. Once the EAG matures it becomes unpalatable for cattle and can cause abscesses in their jaws (Morrow and Stahlman 1984). Identifying the ideal time for grazing in different areas can help ranchers and land managers strategize when it is beneficial to graze, and when it is prudent to avoid grazing a pasture. Grazing practices can help to mitigate wildfire ignitions and spread in EAG invaded areas (Davies et al. 2022(Davies et al. , 2021. Many studies have investigated EAG phenology but only focused on smaller extents (D'Antonio and Vitousek 1992; Davies et al. 2022) like local areas, state or multistate regions, or nature preserves (Bart et al. 2017;Bradley and Mustard 2008;Parker and Schimel 2011;Taylor et al. 2021;Weisberg et al. 2021). These studies have contributed to the knowledge of EAG phenology but lack the larger scale of the western U.S. rangelands. The current study extends knowledge across the western U.S. rangelands that are invaded with EAG and fill the gap for an automated process to help manage and predict EAG growing season for improved management. The study's first objective was to develop a phenology metrics regression-tree model and apply the data from that model to describe the EAG active growing season. The second objective was to increase the understanding of EAG phenology and its relation to weather patterns by analyzing results of a multivariate linear model.

Study area
This research focused on western U.S. rangelands in nine level II Commission for Environmental Cooperation (CEC) ecoregions ( Fig. 1) (Claverie et al. 2018;Dahal et al. 2022). The study area was confined to grassland and shrub land cover defined by the National Land Cover Database (NLCD) 2016 (at the time of processing, the 2019 products were not released), and at or below 2250-m elevation to account for a lack of model training points above that threshold ). Our study area includes more than 66% of conterminous U.S. shrub and herbaceous land cover. For a more detailed investigation, we scale up our results to level III CEC ecoregions to accommodate different ecoregion climates (Wiken et al. 2011).

Data and materials
In order to extract NDVI profiles from likely EAGdominated pixels, we restricted our training points to pixels covered by 20% or 25% EAG abundance using the 30-m Fractional Estimates of Multiple EAG Species datasets (https:// doi. org/ 10. 5066/ P9GC5 JVG) previously predicted by USGS for 2016-2020 (Dahal et al. 2021). Each year of fractional estimates for EAG species was used to determine the potential pixels from which to extract an EAG NDVI profile for manual interpretation. The HLS 30-m NDVI datasets used to extract EAG NDVI profiles included weekly composites from 2016 to 2021 Pastick et al. 2020). Each year consisted of a 52-week NDVI time series profile of a pixel with the requisite EAG estimated fractional cover. We selected training points with over 20 or 25% EAG cover (depending on the HLS tile) because these areas are suitable for EAG phenology detection with a high AMP and distinct differences between senesced and green leaves (Taylor et al. 2021).
The HLS NDVI composites used for extracting NDVI profiles were developed with predicted NDVI values that were compared to actual NDVI values with a strong coefficient correlation (r = 0.79-0.95 for 2017-2021). Fewer HLS satellite scenes in 2016 (about 50% less than each year from 2017 to 2021) reduced the correlation to 0.47  2026 T. D. Benedict et al. (for more information on the development of HLS NDVI composites, see Dahal et al. 2022). While inspecting the HLS NDVI profiles, if an unexpected profile occurred (meaning the profile did not display a usual EAG NDVI profile , also see section NDVI threshold-based method), Google Earth imagery was used for visual interpretation of that pixel's coordinate location to better ascertain the likely vegetation type (we removed the profile if the pixel did not appear to be EAG). Once the manual interpretation of NDVI profiles was complete, our phenology training datasets were used to develop a model for western rangeland EAG phenology metrics. The model also ingested the annual HLS NDVI weekly composites, digital elevation model (DEM) data from the National Elevation Dataset (Gesch et al. 2018), aspect and slope data derived from the DEM, and potential annual incident direct radiation (PADR) (McCune and Keon 2002).
The phenology metrics developed for western U.S. rangelands were then correlated with Parameter-elevation Regressions on Independent Slopes Model (PRISM) (PRISM Climate Group 2021) weather data. The raw weather data had a spatial resolution of 800-m for monthly maximum and minimum temperature (Celsius [°C]) and precipitation (millimeters [mm]) from 2016 to 2021. The 800-m spatial resolution was downscaled to 375-m using the cubic convolution resampling technique to extrapolate details for comparing with spatially aggregated phenology metrics (30-375-m). The aggregated 30-m phenology metrics were resampled by first creating a spatially averaged value within an 11 × 11 (330-m × 330-m) neighborhood around each pixel and then resampling that pixel to 375-m resolution via nearest neighbor. (HLS) tiles separated by extracting training pixels with exotic annual grass (EAG) fractional cover of at least 20% (red) or 25% (blue) Extracting exotic annual grass phenology and climate relations in western U.S. rangeland…

NDVI threshold-based method
To identify EAG sustained growing season's phenological dynamics, we first declared our NDVI thresholds used for both start of season NDVI (SOSN) and end of season NDVI (EOSN). These thresholds served as the foundation for manually extracting phenology metrics based on evaluated NDVI profiles, which were used to build a training dataset of phenology metrics in EAG invaded areas. The training dataset was incorporated into a machine-learning framework, along with other environmental drivers, to extrapolate phenology dynamics across the entire study area.
Developing NDVI thresholds for eight (one of the original nine ecoregions (Marine West Coast Forest) had insufficient NDVI data for threshold extraction) level II CEC ecoregions ( Fig. 1) (Wiken et al. 2011) brought various challenges. Different precipitation, temperature, and biodiversity created environments with diverse EAG growth characteristics in each ecoregion (Prevéy et al. 2014;Thill et al. 1984;Weisberg et al. 2021). Generalizing the start of sustained growth NDVI for EAG in each ecoregion was important to capture various NDVI profiles and their latitudinal distributions. Our NDVI profiles were scaled from the original − 1 to 1 NDVI to 0-200 scaled NDVI using Eq. 1. Table 1 identifies eight ecoregion-specific NDVI thresholds used to determine when (what week) the NDVI profile intercepts the threshold for start of season time (SOST) and end of season time (EOST) (the NDVI is distinctively higher in the growing season (1) Scaled NDVI = (NDVI × 100) + 100 than the dormant season (Taylor et al. 2021)). The NDVI thresholds for SOSN and EOSN (Table 1) were determined by sorting the first set of 7066 randomly stratified training points by calendar year (2016-2020) and their respective level II CEC ecoregion. The stratified training points were used to extract NDVI profiles to manually filter outliers from expected NDVI profiles for EAG. Several criteria were used to determine if an NDVI profile represented an expected EAG . SOST occurred when NDVI had four subsequent weeks of consecutive increasing NDVI values (Taylor et al. 2021), and NDVI values for the start of the five consecutive weeks were greater than or equal to 120 (unscaled NDVI of 0.2) to represent grassland NDVI (Brown 2018). NDVI profiles tended to decrease prior to the first week of July, which, for most if not all the study area, would have been out of the growing season for EAG Prevéy et al. 2014). These criteria capture the NDVI patterns of EAG during SOST and its growing season Huete and Jackson 1987;Taylor et al. 2021). Outliers were confirmed via Google Earth imagery and removed if they included NDVI profiles associated with agricultural fields, large tree canopies or forested areas, and a visually high shrub to grass ratio within a 30-m × 30-m area of the pixel's coordinates. The temporal separation between years helped distinguish the variable phenological cycles that occur each year, especially in arid and semiarid ecoregions (Weisberg et al. 2021;Xie et al. 2022), and can help identify EAG growth patterns that were used to develop SOSN and EOSN thresholds.
The residual NDVI profiles from each calendar year were used to estimate weekly time series of sustained growth. The ecoregion median NDVI values were calculated for each week within the green-up time series, and then the median value was calculated Because our threshold results were so strongly associated with these independent methods, we were confident using our NDVI thresholds to calculate threshold-based phenology metrics.

Manual phenological metrics
The foundation for manual interpreted EAG phenology metrics depended on NDVI threshold values (Table 1) for identifying SOST. Additionally, a list of rules for identifying phenological training datasets was developed to manage similar manual decision making from ~ 11,000 random NDVI profiles (phenology metrics are defined in Table 2). In addition to the rules defined in the previous section, if duplicated values were within the first five weeks of start of sustained growth, that intersected week was rejected as the SOST and subsequent weeks were evaluated for potential SOST. An exception to this rule was when values from the fourth and fifth weeks were equal; if true, then that point would be accepted if the sixth week's NDVI value was greater. An additional exception was applied to pixels failing the previous rules, but the fourth week was greater than week five by at most two NDVI values (or 0.02 unscaled NDVI). If the NDVI profile associated with a training point resulted in a rejected SOST, and there was no other potential SOST, then that point was omitted from the training data. Once the SOST was established for a training point, the senesced period was identified by examining the pattern of descending NDVI values prior to intersecting the NDVI threshold. The EOSN value was identified when the NDVI value met the threshold, or the last week of a consecutive descending pattern, which also determined the EOST. Furthermore, in the senesced period, if successive identical values were at the end of the descending values, the first of the successive values was chosen as EOST because this would represent the end of the sustained growth. In cases of descending then ascending values, we took the lowest value before the second ascending series unless it occured prior to the beginning of July.
Once the SOST and EOST were established, the maximum time (MAXT), maximum NDVI (MAXN), and duration (DUR) were calculated. If the MAXN value occurred more than once in the period between the SOST and EOST, then the last occurrence was recorded as the MAXT. The DUR was determined by the number of weeks between SOST and EOST. Using the above guidelines based on SOST, 36% of the 11,000 NDVI profiles, or 4014 points, were accepted to build the training dataset of seven phenology metrics (SOST, SOSN, EOST, EOSN, MAXT, MAXN, and DUR) for modeling EAG sustained growth.

Modeling and mapping
The EAG phenology metrics training points were extrapolated to the entire study area using predictor  2002). Together, the dynamic and static predictor variables created a model database.
The model database was randomly divided into five subsets with each subset used for validating a model once, while the remaining four subsets were used for model training. The process was repeated four times to build an ensemble of five regression-tree models. This ensemble was developed using Python scikit-learn (Pedregosa et al. 2011) and XGBoost software libraries (Chen and Guestrin 2016). A hyperparameter optimization method, GridSearchCV, was used with an early stop approach to prevent overfitting and to identify optimal parameters during each model iteration. With scikit-learn, we utilized multioutput regressor functionality to estimate six individual phenology metrics (i.e., SOST, SOSN, MAXT, MAXN, EOST, EOSN). We also calculated the number of weeks between estimated SOST and EOST for the DUR metric, and the difference between estimated SOSN and MAXN for the AMP metric. Each set of models was validated using corresponding test data, and results were combined before reporting. Each calibrated model was applied to the spatial versions of the predictor datasets to develop 30-m maps of each metric, making five maps for each metric. Median values and variances of each metrics' five maps served as the final phenology metrics maps and confidence level maps of each metric for 2017-2021. Data for 2017-2020 were developed using manually interpreted training points from their respective years, but 2021 was developed from unseen NDVI datasets to test the robustness of the phenology models.

Weather interaction
We compared EAG species' phenology to monthly weather variables (minimum temperature, maximum temperature, and accumulated precipitation) that may be related to phenology patterns. Precipitation was measured as the amount of rain and melted snow in the given month (PRISM Climate Group 2021) and measured accumulatively from fall to spring (October through December and October through March). The 375-m pixels used for analysis were selected based on their resampled EAG cover of greater than 20% for each year (2017-2021). The EAG temporal metrics (SOST and EOST) were separated based on level III ecoregions' accumulated precipitation prior to the month for SOST starting in October. If the SOST fell on or before the second week of a month, we used the accumulated precipitation from the previous month. The threshold for separation was based on averaging the five-year median accumulated precipitation for each ecoregion, which resulted in two groups: one with at most 200 mm (arid) and the other with greater than 200 mm (semiarid) of precipitation. Ecoregions within level II labeled as Warm Desert (10.2), Cold Desert (10.1), South (9.4) and West Central Semiarid Prairies (9.3) were in the arid group; and the remaining ecoregions from Table 1 were part of the semiarid group. Furthermore, a multivariate linear model was produced (in the statistical software R (R Core Team 2022)) using monthly minimum temperatures, maximum temperatures, and accumulated precipitation from October through December and October through March. When calculating the multivariate linear model, the same variables were used for EOST as SOST with the addition of June's temperature (one of the latest occurrences modeled for EOST). The purpose for these models were not to reproduce a model to estimate phenology metrics, but to identify relations and influences between monthly weather interactions and the SOST and EOST for EAG characteristics.

Phenology modeling results
The EAG phenology model describes 16 EAG species' phenology characteristics. The five iterations of 20% randomly extracted samples used for cross validation verified the robustness of the phenology models. The results from the test points for cross validation showed that the MAXN had the highest correlation (r) of 0.95 followed by the MAXT with an r = 0.94. EOSN and SOST had the next strongest correlations, both with r = 0.92. Figure 2 identifies the model's test scatter plots with the median accuracy between the predicted and observed values for SOST, SOSN, EOST, EOSN, MAXT, and MAXN. Further regression statistics associated with Fig. 2 can be found in Online Resource 2 Table S1. These scatterplots for SOSN and EOSN show the limits the NDVI threshold placed on these phenology metrics based on the Table 1 threshold with the lowest threshold of 126. Figure 3 partial dependence plots (PDP) summarize the top three (out of 42 predictor variables) most important independent variables utilized to develop the phenology model. These plots also show the model's dependencies from each iteration for each phenology metric. The model used week 13 (last week of March) NDVI values most frequently (32% of the time) when determining the SOST for a given pixel, followed by week 14 (20%) and week 12 (11%). Although the model utilized week 13 NDVI as a predictor variable, the model's five iterations estimated areas with lower NDVI in week 13 that were likely to have a later SOST. The five iterations had similar  SOST estimates when the NDVI values in weeks 12, 13, and 14 were between 120 and 140, but fluctuated more when NDVI values were greater than 140. For EOST, the independent variable most frequently used was week 23 NDVI values (32% of the time). The model utilized week 23 and 24 NDVI values to help identify whether EOST will reach the week 27 threshold. Results showed that if the NDVI in week 23 was greater than 135, it was likely that EOST would be week 27. The EOST week 27 threshold is also shown in the EOSN PDP where the top three independent variables were NDVI values from week 26-28, indicating EOSN can be estimated around week 27 values. MAXT had the second strongest correlation (r = 0.94), and estimates were most strongly influenced by weeks 12 (model usage = 20%) and 11 (model usage = 18%). The likely estimate for MAXT occurrence was between week 20 and 21, which was estimated by utilizing NDVI values from weeks 11-13 in all five iterations for modelling MAXT.

EAG phenology metrics
Phenology metric model algorithms were developed and applied to mapping software to generate spatially explicit estimates of six phenology metrics for 16 EAGs (https:// doi. org/ 10. 5066/ P93M8 TEK; Benedict et al. 2022). These phenology metrics are trained based on EAG dominated pixels but are mapped across all vegetation types found throughout the region (Online Resource 1 Fig. S1-S7). Regions of high confidence for SOST occurred within areas that were estimated to have higher fractional cover of EAG. One example of this area is in the Central Basin and Range ecoregion (Fig. 4), where most of the pixels have low EAG cover except for the north-central part (greater than 25% EAG cover) of the ecoregion. SOST estimates in this area were earlier in lower elevation (five-year median SOST around weeks 9-13) than most of the region (week 15). Because the area with a later SOST (week 15) within the Central Basin and Range ecoregion has lower EAG cover (majority of the pixels had 1% cover in 2020), those values could be predicting later native vegetation SOST.
The earliest values for SOST (Fig. 5) that were modeled for the 16 EAG species occurred in the last week of December (week 0) of the previous year. These earlier sustained growths were predicted in the Central California Valley with high confidence. Throughout all ecoregions, the range of median values for 2017 SOST ranged from week 3 (Central California Valley) to week 17 (High Plains The EAG fractional cover shows the estimated percent cover for 2020. c The confidence map shows the pixel-by-pixel confidence for the 2020 SOST phenology metric and Southwestern Tablelands). The median values occurred later in 2018 SOST where the range was week 6 (California Valley) to week 18 (High Plains). These two ecoregions continued to contain the minimum and maximum SOST for 2019, 2020, and 2021 (California Valley SOST week 4, week 3, and week 6 to High Plains SOST week 18, week 17, and week 18, respectively). The median SOST was week 15 (~ second week of April) for the entire study area, but the range of medians throughout the study area was from week 3-18. The median results for 2021 ecoregions were similar to 2017-2020, with values for SOST ranging from week 6 to week 18, which coincided with the other four years of metrics.
The median MAXT for all ecoregions was week 22. The earliest MAXT was week 14 in the Central California Valley and the latest MAXT was week 24 in the Canadian Rockies ecoregion for all five years (Online Resource 1 Fig. S1). In some areas, the MAXT can be the same as the EOST; in these cases, the sustained growth of the 16 EAG species may have The earliest median EOST occurred in the California Valley ecoregion for all five years (2017-2021) at weeks 21, 21, 20, 20, and 19, respectively (Online Resource 1 Fig. S2). The week 27 cutoff for EOST was met in Northwestern Glaciated Plains for all five years. Most of the study area pixels in this ecoregion had 4% cover of EAG in 2020, which would model a later EOST because native perennials are continuing to grow. Only two other ecoregions met the week 27 EOST threshold for four out of five years, the Northwestern Great Plains and High Plains, with most of the study area pixels having 7% and 1% fractional cover, respectively. Overall, the median EOST for each year for all level III ecoregions occurred between weeks 24 and 25 (Online Resource 1 Fig. S2).

Weather interaction results
The multivariate linear relations between temperature, accumulated precipitation, and EAG phenology are complex and variable, especially across broad geographic areas. Both negative and positive relations occur, many times with the same variable from one year to the next. Weather patterns in arid and semiarid environments can be highly variable annually, so the change in relations between weather and phenology is expected. From an ecological perspective, patterns should emerge and generally be comprehensible when comparing relations between weather and EAG phenology. One note is that the years, e.g., 2017, represent the data for the growing season, e.g., 2016-2017, so the October (suffix 10) in the variable column represents weather from October 2016.
In ecoregions with the average precipitation prior to SOST greater than 200 mm (Table 3), the strongest relation for SOST was negative, when the lower 2017 May temperature maximum induced a later SOST for EAGs in that growing season. Other strong negative relations included accumulated precipitation from October 2018 through December 2018, which indicated that the more precipitation received, the earlier SOST occurred in the 2019 growing season, and the lower the 2018 and the 2021 growing seasons' October temperature maximums, the later the SOST occurred in those two growing seasons. The strongest of the positive relations was temperature maximum for January 2020, which indicated that warmer temperatures led to later SOST. October temperature maximums for both the 2017 and 2020 growing seasons had two of the weakest relations with SOST compared to the remaining years, signaling other variables (i.e., 2017 May temperature maximum and 2020 January temperature maximum) were more influential on their respective years' SOST. Because many EAG species are winter annuals that, with adequate precipitation, can germinate and start growing in the fall, October can be a crucial month for EAG SOST in these ecoregions.
Because most of the variables between temperature and accumulated precipitation for semiarid ecoregions showed more significant relations (t-statistic was greater than one standard deviation from the mean) with maximum temperatures, ecoregions with less than 200 mm (arid ecoregions) of accumulated precipitation prior to SOST (Table 4) had minimum temperatures play more of a role in its significant effect on SOST. Of all temperature variables, December 2016 temperature minimum affected SOST for the 2017 growing season with a strong negative relation, which indicated warmer minimum temperature brought on earlier SOST. Negative relations held true for all years' SOST and minimum temperatures in December except for the 2020 growing season. The strongest negative influence was in the 2021 growing season between accumulated precipitation from October through December 2020, which indicated more precipitation in fall, the earlier SOST. Both December's minimum temperature and accumulated precipitation showed the same relations to SOST with the exception of 2018. In the 2018 growing season the most significant positive relation occurred with temperature maximum in December 2017. Other than the 2018 multivariate linear model, the other years' SOST most influential variable (t-statistic was greater than one standard deviation from the mean) occurred during the winter months.
The SOST is one of the more important phenology metrics, but another important metric describing EAG growing season is the EOST. The arid ecoregions (Online Resource 2 Table S2) were negatively affected by the minimum temperature in March in all five years for EOST. When March temperatures were strongly related to EOST, the maximum and minimum temperatures also affected and decreased the precipitation effectiveness for predicting EOST for 2018 and 2019 growing seasons. The minimum Table 3 Five multivariate linear models for ecoregions with greater than 200 mm average accumulated precipitation prior to exotic annual grass (EAG) start of season time (SOST) from 2017 to 2021 The variables include maximum temperature (tmax), minimum temperature (tmin), and accumulated precipitation (Accppt) starting from October (10) Table 4 Five multivariate linear models for ecoregions with less than 200 mm average accumulated precipitation prior to exotic annual grass (EAG) start of season time (SOST) from 2017 to 2021 The variables include maximum temperature (tmax), minimum temperature (tmin), and accumulated precipitation (Accppt) starting from October (10) from the previous year, e.g. Accppt10-03 2018 is precipitation from October (10) 2017 through March (03) 2018. The "estimate" columns are the coefficients calculated for the respective variable (temperature in Celsius and precipitation in millimeters) followed by the standard error (std.error), t-value from the t-test (statistic), p.value where p-values < 2E-16 are represented as "***", and the coefficient of determination (R 2 ) for each model.  temperatures in March had more of an effect on EOST than the maximum temperatures based on its coefficients. If higher minimum temperatures in March occurred, and the maximum temperature remained consistent, the linear models indicated that EAG species had an earlier EOST. For both semiarid and arid ecoregions (Online Resource 2 Table S3 and S2, respectively), during the 2017 through 2019 growing seasons, EOST was more influenced by temperature than accumulated precipitation. For the 2020 and 2021 growing seasons, temperature had more effect on EOST, but the relations compared to accumulated precipitation showed a strong negative influence between the accumulated precipitation through December and a positive relation through March.

EAG phenology model
We developed an automated process to estimate EAG sustained growth characteristics by creating a model from NDVI threshold-based phenology training data. We considered that EAG species have variable phenological patterns year to year (Bart et al. 2017;Bradley and Mustard 2008) and manually identified the sustained growing season for each year using annual NDVI phenology profiles. Grouping 16 EAG species' phenology characteristics can capture alternating growing seasons, as mentioned in Weisberg et al. (2021). For example, cheatgrass in its senescent phase, can occur concurrently with medusahead in its active growing stage. The strong correlations (r = 0.88-0.95, Fig. 2) from the five iterations of cross-validation indicate robust multi-output models for estimating all six EAG phenology metrics on unseen data and through time .
Modeling natural occurrences using remotely sensed data and machine learning is challenging and comes with uncertainties. One challenge we identified in this study was creating a one-size-fits-all metric and model throughout several ecoregions. For example, EAG pixels in the level II ecoregion Mediterranean California (Wiken et al. 2011), showed HLS NDVI value differences greater than 20 (unscaled NDVI difference of 0.2) for week 52 of one year and week 1 of the next. The large difference implies a substantial increase in greenness in this one-week period, which does not seem likely for a natural environment and may be affected by predicted HLS data for nonexistent NDVI composites . However, some studies have recorded cheatgrass SOST beginning between November and December Schimel 2010, 2011). Also, Bolton et al. (2020) showed 50% green up in areas of California starting prior to January 14, indicating earlier SOST for EAG near the Central California Valley.

EAG phenology metrics
In all five years of the study period, areas with less than 10% EAG cover had later SOST than areas dominated by EAG. Bradley, Mustard (2008) found the mean SOST for cheatgrass in the Great Basin area was day of year (DOY) 92 (first week of April) and the mean range was 64-111 (from the first week of March to the third week of April) for 1991-2000. Our results for 16 EAG species (cheatgrass being the most dominant species Dahal et al. 2022;Pastick et al. 2020) in the Great Basin (level III ecoregions Snake River Plain, Central Basin and Range, and Northern Basin and Range) average SOST was week 13.7, which was within a week of Bradley & Mustard's (2008) cheatgrass metrics.
Other phenology metrics from this study compare well with prior work. The median value of MAXT, calculated utilizing all ecoregions, was week 22 (last week of May). This was similar to the Bradley and Mustard (2008) study where they calculated a mean MAXN on DOY 149 (within the last week of May) and their mean range was DOY 126-164. The MAXN range in the 2021 growing season was comparable to other years in the study period and shows the robustness of our phenology model. Pixels unseen by the model algorithms returned similar phenology metrics as previous studies. Understanding the timing of sustained growth and its dynamics for EAG can help land managers understand when to use grazing techniques that can reduce EAG abundance and provide access to early season grazing resources (Davies et al. 2022). Ranchers can target grazing to EAG invaded pastures during the peak growing season (Browning et al. 2019). Thill et al. (1984) found that in the Idaho Batholith, cheatgrass started "heading" from late April to early May. For cheatgrass to head (or flower) in late April indicates the species began growing prior to flowering and therefore the start of season was before week 17 (Thill et al. 1984). Results from our study show that the median SOST in the Idaho Batholith was during weeks 15-16 (early to mid-April) across four years. The general agreement between this phenology study's data and previous independent studies' data provides a level of confidence in our phenology metrics.

EAG weather variables
Many variables can affect the timing of the germination, growth, and senesced stages (D'Antonio and Vitousek 1992) of a species, and examining the multivariate linear models in the level III ecoregions can give us a better understanding of EAG lifecycles. We chose to examine weather data as it relates to EAG species because weather, especially precipitation, is a strong driver of EAG growth dynamics. These multivariate linear models were developed to identify the influence that monthly weather data may have on EAG phenology. We noticed that in wetter ecoregions, SOST was easier to predict using only monthly weather variables (average coefficient of determination (R 2 ) = 0.55). These multivariate linear models treat the variables as independent from each other, but overall influence the dependent variable as a whole. These influences consider the differences that each variable has on the dependent variable (SOST or EOST) and identifies a significant relation, either negative or positive. The variability of these relations likely occur because EAG life cycles are dependent on multiple driving forces, including less dominant temperature and precipitation drivers, as well as soil nitrogen (D'Antonio and Vitousek 1992; Parker and Schimel 2011), soil moisture (Bart et al. 2017;D'Antonio and Vitousek 1992;Thill et al. 1984), and seed depth (Thill et al. 1984).
Our data corresponded well with the Idaho Batholith (one of the wetter ecoregions) results mentioned in Thill et al. (1984). They observed winter-dormant cheatgrass restarting growth after being dormant when temperatures warmed, likely capitalizing on soil moisture that was a product of winter and/or early spring precipitation. The arid ecoregions analyzed showed strong correlations between accumulated precipitation and SOST and EOST, which agrees with prior research conducted in Cold Desert ecoregion and other arid regions (Weisberg et al. 2021;Xie et al. 2022). In the Mediterranean California ecoregion (a wetter ecoregion), Bart et al. (2017) found that EAG phenology was moisture dependent. They also found that the senesced period of Mediterranean grasses in two study sites (both in the California Coastal Sage ecoregion) was dependent on multiple variables, including temperature and the timing and amount of precipitation (Bart et al. 2017). Similar patterns were found in semiarid ecoregions' linear models, but for the 2020 growing season SOST had less significant relations to precipitation. In contrast, the 2020 growing season EOST showed more of a significant relation with precipitation.
The EAG metric models identified the detectable sustained growing season using NDVI values. Weather patterns and events from periods' preceding the actual onset of EAG germination can influence germination onset and even post-germination growth, with freezing temperatures stopping EAG growth and growth progressing when temperatures warm (Bart et al. 2017;Thill et al. 1984). Other driving effects for germination besides weather variables could be where monthly weather had a less significant effect on SOST or EOST. Cheatgrass is one of the dominant species within the 16 EAG species modeled but is not the only one being measured; thus, the extended EOST may take into account the other species' growing season showing a longer EOST and other weather affecting the end of season. Because cheatgrass is one of the most studied EAG species, and other EAG species may respond differently to monthly weather estimates than cheatgrass, these linear models may indicate different significant variables than found in other studies. Weisberg et al. (2021) observed these overlapping growing seasons with cheatgrass and medusahead; therefore, gaining more precipitation after December may prolong the EOST and not end as soon as one may expect with cheatgrass.

Conclusion
The main foci of this study are twofold: (1) to describe the EAG active growing season and (2) to generate a more in-depth understanding of the relations between weather patterns and EAG phenology. The foci help to identify when sustained growth of EAG begins and ends in western rangelands and to understand the phenology dynamics in this highly variable landscape. The study's results can be used by researchers, land managers, and conservation teams to develop plans that benefit their needs for controlling or utilizing (grazing) EAG species. Utilizing the temporal metrics such as SOST, MAXT, or EOST can benefit land managers by identifying preparation periods for wildfire danger, cattle grazing, or herbicide applications. With EAG spread throughout rangeland ecoregions, identifying when these species are flammable (EOST) can help identify vulnerable locations that may benefit from preventive measures. If a grazing technique is utilized, SOST can also help land managers to identify when to move the livestock into rangeland areas invaded by EAG species for grazing prior to seed production to ensure the best nutrient rich product for cattle. Knowing EOST or MAXT can help ranchers time the removal of livestock prior to seed maturity to prevent EAG spatial distribution and harmful seed ingestion from some EAG species. Some mitigation plans that use herbicides recommend applying chemicals prior to EAG emergence, so utilizing the SOST time series can identify patterns for crucial application periods. NDVI metrics may be an asset with other remotely sensed data for monitoring areas out of range or prevent foot traffic through invaded areas.
For semiarid ecoregions, we were able to identify the importance that October has on SOST. The paper investigated how these sustained growth characteristics relate to monthly temperature and accumulated precipitation per growing seasons and the difference in variability between semiarid and arid ecoregions. Within these variabilities, some trends emerged (both negative and positive) based on the monthly weather variables throughout the growing season. These patterns were identified within a multivariate model to demonstrate how the months' weather may affect EAG SOST or EOST. Utilizing monthly weather data and their relations from the linear models show their robust responses to these weather patterns. Comparing collected weather data from a current growing season to the phenology metrics may allow a rough prediction on earlier SOST or EOST. Future research may include processing more EAG NDVI profiles in the Mediterranean California (level II) ecoregion to identify earlier SOST and utilizing and fine-tuning our phenology metrics for a more robust model. The EAG phenology metrics time series can increase the understanding of EAG dynamics and help land managers develop plans to protect entities deemed valuable to society like native habitat, wildlife, recreational areas, and air quality.

Conflict of interest
The authors declare no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.