Differentiating Cheatgrass and Medusahead Phenological Characteristics in Western United States Rangelands

Abstract

Expansions in the extent and infestation levels of exotic annual grass (EAG) within the rangelands of the western United States are well documented. Land managers are tasked with developing plans to limit EAG spread and prevent irreversible ecosystem deterioration. The most common EAG species and the subject of extensive study is Bromus tectorum (cheatgrass). Cheatgrass has spread rapidly in western rangelands since its initial invasion more than 100 years ago. Another concerning aggressive EAG, Taeniatherum caput-medusae (medusahead), is also commonly found in some of these areas. To control the spread of EAGs, researchers have investigated applying several control methods during different developmental stages of cheatgrass and medusahead. These control strategies require accurate maps of the timing and spatial patterns of the developmental stages to apply mitigation strategies in the correct areas at the right time. In this study, we developed annual phenological datasets for cheatgrass and medusahead with two objectives. The first objective was to determine if cheatgrass and medusahead can be differentiated at 30 m resolution using their phenological differences. The second objective was to establish an annual phenology metric regression tree model used to map the growing seasons of cheatgrass and medusahead. Harmonized Landsat and Sentinel-2 (HLS)-derived predicted weekly cloud-free 30 m normalized difference vegetation index (NDVI) images were used to develop these metric maps. The result of this effort was maps that identify the start and end of sustained growing season time for cheatgrass and medusahead at 30 m for the Snake River Plain and Northern Basin and Range ecoregions. These phenological datasets also identify the start and end-of-season NDVI values, along with maximum NDVI throughout the study period. These metrics may be utilized to characterize annual growth patterns for cheatgrass and medusahead. This approach can be utilized to plan time-sensitive control measures such as herbicide applications or cattle grazing.

1. Introduction

Expansions in the extent and infestation levels of exotic annual grass (EAG) within the rangelands of the western United States are well documented. The most common EAG species and the subject of extensive study is Bromus tectorum (cheatgrass) [1,2,3,4,5,6,7,8,9,10]. Cheatgrass has spread rapidly in western rangelands since its initial invasion more than 100 years ago [11]. Another aggressive EAG, Taeniatherum caput-medusae (medusahead), is also commonly found in the western United States [12]. It has been known for decades that medusahead can replace cheatgrass and other annual grasses and is an equal threat in those areas with cheatgrass [13]. In addition, EAG species such as cheatgrass and medusahead often have interspecific competitions in western U.S. rangelands. These EAG species have been known to outcompete native shrubs (e.g., sagebrush (Artemisia spp.)) and grasses (e.g., bluebunch wheatgrass (Pseudoroegneria spicata) and Sandberg bluegrass (Poa secunda)) and even displace wildlife (e.g., the greater sage-grouse (Centrocercus urophasianus)) in sagebrush ecosystems [7,14,15]. Land managers are tasked with developing strategy plans to slow EAG spread and prevent irreversible ecosystem deterioration. Many studies have focused their efforts on mapping cheatgrass [3,6,9,16,17,18] and medusahead distributions [12,17,18,19,20], providing important spatial information for management. In addition, the phenological patterns associated with these EAG-occupied areas are important for managing invaded ecosystems.

For EAG, phenology patterns can change year to year depending on seasonal temperatures and precipitation [21]. Identifying when these changes occur for a location is important for monitoring changes in biodiversity or degradation, and remote sensing technology makes monitoring larger areas more feasible when considering costs and effort. Areas dominated by native perennial species generally have longer growing seasons than areas dominated by EAG species, and a change in phenological patterns after a disturbance may signify the beginning of an invasion or change in the environment [22,23]. Applestein and Germino [24] noted a positive correlation with cheatgrass and medusahead invasions after fires in sagebrush ecosystems. The longer the period after a fire disturbance, the higher the probability of cheatgrass invading followed by medusahead. Early identification of these disturbed areas allows time for treatments to be planned and implemented to prevent cheatgrass and medusahead from dominating the land cover. A better understanding of EAG phenological patterns can help land managers plan effective treatments to limit the spread of these invasive species.

Phenological timing is important for the management of these short-growing-season grass species and can make a difference between unsuccessful and successful efforts [25,26]. Some management practices require time-sensitive application for successful eradication or control. Rinella et al. [26] studied the effectiveness of herbicides for managing cheatgrass and medusahead during distinct plant growth stages in California. Herbicide applications were effective when the herbicide was applied during medusahead’s preheading stage, which was also noted as the easiest time to differentiate medusahead from cheatgrass [26]. Other management actions included mowing or grazing to control spreading. Brownsey et al. [25] studied the effects that mowing and grazing at different phenology stages had on medusahead control. The results showed that grazing was only effective during the livestock palatability stage defined as the “spike emergence, crude protein, and acid detergent fiber” in the plant [25]. Furthermore, they identified the boot stage, which normally lasts for 10 to 15 days, as the time in which EAG contains nutrients suitable for grazing and impacts seed development [25]. Phenology products could inform management decisions to use livestock for medusahead control during the boot stage.

As is common with winter annual grass species, the timing of the phenological stages of cheatgrass and medusahead tends to vary from year to year. The difficulty in mapping these species separately is due to their co-occurrence in moderate spatial resolution imagery, similar spectral profiles, and distribution. Weisberg et al. [18] used an unoccupied aerial vehicle (UAV) to map 2.78 ha intermixed with cheatgrass and medusahead. One of their goals was to determine the spectral difference in cheatgrass and medusahead utilizing visible reflectance and near-infrared. They found the visible reflectance was valuable in discriminating the two species throughout their growing season, which provided accurate classifications at a fine-scale resolution (2 cm). Weisberg et al. [18] were able to distinguish cheatgrass and medusahead spectrally at fine resolution, but coarser resolution may be difficult. Clinton et al. [27] were able to use MODIS 250 m normalized difference vegetation index (NDVI) 16-day composites to map cheatgrass abundance. Although Clinton et al. [27] did not differentiate cheatgrass and medusahead, they showed it was possible to map cheatgrass at a coarser resolution. Because medusahead generally matures later than cheatgrass, it may be possible to identify each species using their phenological characteristics [18,24]. Hironaka [13] noted medusahead matured 2 to 3 weeks after cheatgrass matured in Gem County, Idaho. Capturing the different growing seasons of cheatgrass and medusahead may require moderate- to high-temporal-resolution remotely sensed products. One product that fulfills the temporal resolution is weekly (7-day) NDVI composites derived from Harmonized Landsat Sentinel-2 (HLS) instruments [28]. NDVI is one of the most widely used vegetation indices for phenology based on its availability and ability to monitor greenness and vegetation health [29]. Qin et al. [29] compared leaf area index (LAI) to multiple indices used for extracting phenology; a few of these indices are NDVI, enhanced vegetation index (EVI), and normalized difference phenology index (NDPI). LAI has been closely related to photosynthetic activity, which can be comparable to phenology [29]. NDPI and NDVI were in the top three indices correlated to LAI in both Landsat 8 and Sentinel comparisons, suggesting them to be good indices for extracting phenology. Additionally, HLS weekly NDVI has been proven effective for monitoring the distribution of cheatgrass and medusahead [6,10,17,30]. Therefore, the temporal and spatial resolution of HLS NDVI may differentiate cheatgrass’s and medusahead’s growing season differences.

This work has two research objectives to improve understanding and develop spatially and temporally extensive geospatial products for managing EAG species. The first objective was to determine if extracted phenological signals of cheatgrass and medusahead can be differentiated at 30 m resolution. The second objective was to establish an annual phenology metric regression tree model and map the growing season of cheatgrass and medusahead. This research promotes the innovations of new datasets to help increase the scientific knowledge of the spatial distribution of cheatgrass and medusahead phenological characteristics. Land managers and researchers may use these metrics to identify annual growth patterns of cheatgrass and medusahead. These products can be utilized to plan time-sensitive control measures such as herbicide applications or cattle grazing for maximum nutrition content and EAG control.

2. Materials and Methods

2.1. Study Area

The Snake River Plain (SRP) and the Northern Basin and Range (NBR) ecoregions [31] were selected as the study area due to the prevalence of cheatgrass and medusahead according to their respective historical fractional cover from U.S. Geological Survey (USGS) estimates [30]. The SRP and NBR (Figure 1) are located in the level II Commission for Environmental Cooperation (CEC) Cold Desert ecoregion based on their environmental characteristics [31]. From Parameter-elevation Regressions on Independent Slopes Model (PRISM) data for 2016 through 2021, SRP and NBR experienced a monthly average of 34 mm and 37 mm precipitation during the germination months (October through March), respectively [32]. Temperature for germination months ranges from an average minimum of −3.9 Celsius (°C) to an average maximum of 7.2 °C [32]. There was an average of 3 million hectares of cheatgrass and 12 thousand hectares of medusahead in the SRP and NBR based on cover estimates between 2017 and 2021 based on pixels estimated to contain at least 20% cover [20].

Elevation in the study area ranges from 630 m to 3303 m above sea level based on the digital elevation model (DEM) [33]; however, we masked out elevations above 2350 m because of a lack of model training data above this threshold (Figure 1). Furthermore, we masked areas not classified as shrub or grassland by the 2019 National Land Cover Database (NLCD) [9,17,21,34]. These restrictions provide sufficient indications of cheatgrass and medusahead phenology derived from NDVI values.

2.2. Data

Identifying phenology characteristics from an NDVI time series requires appropriate spatiotemporal resolution to identify sustained growth. HLS-based cloud-free weekly 30 m NDVI datasets were used to develop training data of cheatgrass and medusahead sustained growth [17]. For more information on the development of HLS-based weekly cloud-free NDVI composites, refer to Dahal et al. [17]. More than 298,000 scenes from HLS (v. 2.0) with 2–3-day temporal resolution were used to create the weekly (7-day) cloud-free composites for 2016–2021. Utilization of machine learning enabled the development of quality cloud-free NDVI composites with a strong correlation coefficient between predicted and measured NDVI (r = 0.79–0.95 for 2017–2021, r = 0.47 for 2016 due to fewer HLS scenes) [17]. Weekly NDVI composites at 30 m resolution were used in this study because they provide sufficiently fine temporal frequency and cloud-free imagery to study phenological differences in short-lived cheatgrass and medusahead grass species. The weekly composite values formed a continuous 52-week NDVI profile to describe the growth of cheatgrass and medusahead over an annual phenological cycle (Figure 2). The phenology metrics in Figure 2 describe the detectable greenness based on NDVI time series. These metrics describe when cheatgrass or medusahead green up, how long they stay green, and when maximum greenness occurs. The resulting dynamic variables used in the model were 33 weeks of NDVI values (i.e., week 48 NDVI values from the previous year through week 28 NDVI of the targeted year). These weeks were determined to capture the intended growing season time frame for cheatgrass and medusahead [21].

Static variables that were used as inputs in the cheatgrass and medusahead phenology model included digital elevation model (DEM) data from the National Elevation Dataset [33] and derivative products of slope and aspect. Additionally, potential annual incident direct radiation (PADR) was also used as a variable to develop the models [35]. Another set of static variables included the properties of the topsoil, which includes soil organic matter, available water capacity, clay, silt, and sand content [17,36]. Finally, we included 30-year climate normals of annual and winter normals for precipitation, maximum temperature, and minimum temperature from 1988 through 2023 from Daymet [20,37]. The total set of predictor variables included 15 static variables and 33 dynamic NDVI variables.

2.3. Developing Training and Test Data

The methods consisted of a three-step process. We first identified a set of points to derive an HLS NDVI time series based on a high probability of cheatgrass and medusahead cover. Second, we extracted the phenological metrics used for training the models by applying a decision tree processing technique on the NDVI time series. Finally, we utilized automated machine learning techniques to derive phenological models that were used to develop maps for the entire study area per 30 m pixel.

Training and test data for this study were extracted from the HLS NDVI time series and driven by historical annual fractional cover maps of cheatgrass and medusahead that are released via the USGS Rangeland Exotic Plant Monitoring System team [20]. We used a stratified sampling technique to extract random pixels with a high probability of at least 20% individual cover of cheatgrass or medusahead (Figure 3a) [5,27]. To improve the certainty of the high-probability areas, we eliminated clusters of pixels with less than five high-probability pixels within a 3 × 3 cell moving majority window [10] and used the pixels that intersected these two restrictions. The selected pixels were further restricted based on the co-occurrence of cheatgrass and medusahead to create more species-pure training. Specifically, we required pixels to contain at least 50% more relative cover of one species than the other (e.g., a site with 25% cheatgrass and 20% medusahead would not be selected, but a site with 25% cheatgrass and 5% medusahead would) [38]. In 2022, medusahead did not exceed 50% more relative cover than cheatgrass in the same pixels, so random samples were taken within pixels with at least 15% cover instead of 20%. The pixels were further filtered using the confidence layer associated with the cheatgrass and medusahead cover maps (ranging from 0 to 10, with higher values indicating more confidence in the percent cover prediction [17]). Confidence values greater than 0 were required for medusahead and greater than 9 were required for cheatgrass; due to data availability, cheatgrass had a higher confidence restriction. NDVI values were processed to 8-bit unsigned scaled NDVI (scNDVI) values from Equation (1):

$$\mathrm{scNDVI}=(({\rho }_{NIR}-{\rho }_{red})/({\rho }_{NIR}+{\rho }_{red})\times 100) + 100.$$

(1)

The phenology metric processing approach used in this study was developed based on an NDVI threshold approach [39,40,41]. The threshold approach identifies a specific NDVI value that signifies when the growing season occurs within an NDVI time series [42]. The NDVI threshold was calculated twice. Figure 3a shows the first calculation of the threshold used to develop an initial assessment using the stratified samples to determine ideal phenological profiles and eliminate impractical NDVI time series. Figure 3c shows the second calculation for determining the final NDVI threshold used to evaluate the NDVI time series for training models. To capture the NDVI threshold, 10% of the maximum unscaled NDVI value within the growing season (up to week 30) was added to the minimum unscaled value (up to week 24) to determine an acceptable start of sustained increasing NDVI and a sustained decrease in value [39,40,41]. The first calculation for the NDVI threshold was averaged for all random profiles extracted (~150,000 pixels). The final NDVI threshold (126 scNDVI) used to extract training data was developed from accepted profiles after utilizing the first calculated threshold (Figure 3c).

Once the NDVI threshold was determined, it was utilized in the phenology metric analysis. We used decision tree analysis to extract phenology metrics from the ~150,000 stratified random NDVI time series (Figure 3c). First, start-of-season time (SOST) and start-of-season NDVI (SOSN) were found by analyzing when the weekly NDVI time series reached the threshold and if the following five weeks had a consistent increase in NDVI. Next, the end-of-season time (EOST) and end-of-season NDVI (EOSN) were determined based on a steady decrease towards the NDVI threshold. If the NDVI threshold was not reached by week 25, EOST and EOSN were set based on week 25 NDVI values. Week 25 was selected as a training data limit because the climate may be too warm and dry for these species to continue growth past the third week of June [9,43]. Lastly, the maximum time (MAXT) and maximum NDVI (MAXN) were identified within the weekly time series between SOST and EOST. If all values were less than the NDVI threshold, or the profile did not show consistently increasing NDVI values, the pixels were rejected. The accepted profiles from the decision tree analysis were further curated to profiles with at most a MAXN of 172. The MAXN of 172 was chosen based on the maximum MAXN used in the EAG phenology training data to remove pixels that were least likely cheatgrass or medusahead [21]. The remaining accepted profiles were utilized as training data to represent sustained growth for respective high-probability cheatgrass and medusahead cover.

The decision tree rules were tested among the accepted 6591 NDVI profiles that were used to train EAG phenology models from Benedict et al. [21] (Figure 3b). The EAG phenology training data were utilized to test how well the decision tree analysis extracted phenology metrics based on the manually extracted metrics from Benedict et al.’s [21] methodology. We evaluated for a high correlation in Pearson’s r and a low error in mean absolute error (MAE) to verify if the decision tree rules matched the results.

Figure 3. Phenology method flowchart. This flowchart is the overview of capturing the phenology training data for cheatgrass (BRTE [44]) and medusahead (TACA8 [44]). (a) The methods for extracting high probability BRTE and TACA8 pixels, (b) the decision tree analysis development using exotic annual grass (EAG) training data from Benedict et al. [21], (c) developing the training data, and (d) developing the phenology model.

Figure 3. Phenology method flowchart. This flowchart is the overview of capturing the phenology training data for cheatgrass (BRTE [44]) and medusahead (TACA8 [44]). (a) The methods for extracting high probability BRTE and TACA8 pixels, (b) the decision tree analysis development using exotic annual grass (EAG) training data from Benedict et al. [21], (c) developing the training data, and (d) developing the phenology model.

2.4. Modeling Phenology

Utilizing machine learning techniques, Python scikit-learn [45], and XGBoost software (version 2.0.3) libraries [46], we created an ensemble of five regression tree models with a multi-regressor output for SOST, MAXT, and EOST (Figure 3d). XGBoost’s early stop approach based on MAE was used for the calibration of the models. The classification determined pixel-by-pixel values and was not influenced by surrounding pixels. Furthermore, the final estimated maps were the median values of the five regression tree model iterations (Figure 3d). The model was trained on 2017–2021 data and was then applied to 2022 variables and NDVI datasets to develop an annual map that estimates cheatgrass and medusahead phenological patterns.

Using the estimated SOST, MAXT, and EOST, a pixel-drilling method was developed to extract the remaining phenology metrics (i.e., SOSN, EOSN, MAXN, duration (DUR), and amplitude (AMP)) based on the individual pixel’s NDVI raster (Figure 3d). The pixel-drilling method extracted the actual NDVI values per pixel’s NDVI time series at its associated temporal metric for SOSN and EOSN. MAXN was extracted as the relative maximum NDVI between SOST and EOST in the NDVI time series. This was developed to improve the previous EAG phenology model (developed by Benedict et al. [21]) to minimize errors for values that can be extracted from their respective NDVI profiles using the model-estimated temporal metrics.

A five-fold cross-validation approach was applied to calculate accuracy and error within the estimated temporal metrics. Each model was run five times with a unique subset of 20% of the training data withheld as validation, resulting in five cross-validations. Therefore, each data point was used once to test a model and four times to train a model. The final map represents the median values of the five mapped iterations. The final maps were then evaluated using 20% withheld from the original set of training data for testing, data that were not included in model development for cheatgrass or medusahead (Figure 3d labeled as “20% held for testing”). Pearson’s correlation coefficient (r) was used to evaluate how the observed values compared to the estimated values, where a high r value would indicate a highly correlated model. To test the model’s accuracy, we calculated the MAE using Equation (2) and then normalized the root mean square error (RMSE) (Equation (3)):

$$\mathrm{MAE}={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N|y_i-{\hat{y}}_i|,$$

(2)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{(y_i-{\hat{y}}_i)}^2}.$$

(3)

where $y$_irepresents the ith measured value and $\hat{y}$_i represents the ith predicted or estimated value. Lower MAE indicated better model fit. To normalize RMSE (hereafter relative root mean square error (RRMSE) (%)), we divided Equation (3) by the measured average and multiplied by 100 to determine how good of a fit the model was compared to the measured values [47,48]. The standard interpretation was that a model with excellent fit would have an RRMSE of <10%; for a good fit, 10% < RRMSE < 20%; fair, 20% < RRMSE < 30%; and a poor model would have RRMSE > 30% [47,48,49]. The Mann–Whitney U-test was used to test for significant differences between cheatgrass and medusahead due to non-parametric testing and results from the Shapiro–Wilk test showing a non-normal distribution. The rank-biserial correlation was used to calculate the effect size of the Mann–Whitney U-test [50]. Rank-biserial correlations represented proportions of medusahead and cheatgrass, and values ranged from −1 to 1. A value of 0 indicates that cheatgrass and medusahead were present in equal proportions, a value of −1 indicates complete coverage by medusahead, and a value of 1 indicates complete coverage by cheatgrass [50].

3. Results

3.1. Phenology Decision Tree Analysis

The application of the decision tree rules to the EAG phenology training data (6591 NDVI profiles in Figure 3b) resulted in the acceptance of 94% (6191 points) of the training data used for the EAG phenology model from Benedict et al. [51].

Table 1. Pearson’s r correlation and mean absolute error (MAE) for phenology decision tree analysis validation using the 6591 training data points used in building the exotic annual grass phenology model [21]. MAE units for start-of-season time (SOST), end-of-season time (EOST), and maximum time (MAXT) are based on weeks. MAE units for start-of-season NDVI (SOSN), end-of-season NDVI (EOSN), and maximum NDVI (MAXN) are based on scaled NDVI values.

	SOST	SOSN	EOST	EOSN	MAXT	MAXN
Pearson’s r	0.96	0.94	0.89	0.93	0.99	1.00
MAE	0.39	0.65	0.92	1.79	0.08	0.03

shows the correlation and MAE from comparing the phenology decision tree output and the EAG training data. Overall, there was a high correlation for each metric. The EOST metric had the lowest correlation (0.89) and the highest MAE within the temporal metrics, although the MAE was still less than 1 week. The EOSN had the largest error for NDVI metrics (MAE = 1.79 scNDVI values) and the second lowest correlation coefficient (0.93) [21].

3.2. Phenology Training Data

Quality training data are crucial for developing useful models, and understanding the training data is helpful in managing expectations and understanding data characteristics. While evaluating the training data, some relations were identified in the SOST and MAXT between cheatgrass and medusahead. In Figure 4, the boxplots show the distributions of the training data for years 2017–2021 with a median SOST for medusahead equal to week 17 and cheatgrass equal to week 12. The MAXT had a smaller difference as medusahead reached its maximum NDVI three weeks after cheatgrass on average. Although the yearly median SOST changes from year to year for both medusahead and cheatgrass, we observed a phenological lag every year where medusahead SOST occurs after the cheatgrass SOST, consistent with expectations from the literature. In 2021, there was more overlap in the SOST, where the first quartiles of cheatgrass and medusahead were equal. The 2021 medians also were the latest SOST for cheatgrass and earliest SOST for medusahead. Knowing this, we would expect 2021 phenology metrics to begin earlier for medusahead when compared to previous years as seen in Figure 5.

Comparing medians between cheatgrass and medusahead, there seems to be a difference in the extracted NDVI profiles. Both species profiles were subject to the same NDVI threshold and the same phenology assessments, resulting in divergence between species.

3.3. Model Results

The five-fold cross-validation results for both species were calculated using the median between the five modeling iterations. Overall, observed and estimated weeks of medusahead SOST and MAXT were highly correlated and had low MAE, whereas correlations between observed and estimated weeks of medusahead EOST were relatively high but were the lowest of all regression models. As shown in Figure 6, there are some deviations from the 1:1 line, but most closely followed the observed-equals-estimated relation. Cheatgrass models used 31,281 data points, and models for medusahead included 14,798 data points, which may have contributed to the lower Pearson’s r values for the EOST model for medusahead. The SOST and EOST values for the medusahead (Figure 6d,e) model occurred later in the season and had a smaller range compared to SOST and EOST values for cheatgrass (Figure 6a,b). The medusahead model showed a better accuracy than the cheatgrass model with MAE less than a quarter of a week off for temporal metrics, whereas cheatgrass MAXT and EOST MAE were more than a third of a week off.

The cross-validation statistics are important for estimating temporal metrics, but the maps must match expected patterns to be useful. We inspected the model-estimated outputs using the residual 20% phenology profiles that were not used in the model training (Figure 3d “20% held for test”). These points validated cheatgrass and medusahead maps for areas with at least 20% cheatgrass and 20% medusahead cover, respectively. Pixels with at least 20% cover of cheatgrass showed a high correlation (r > 0.85) for each year’s phenological metrics when comparing the phenology decision tree results with the modeled results (

Table 2. Cheatgrass (BRTE [44]) and medusahead (TACA8 [44]) estimated phenology metrics compared to decision-tree-produced phenology. The estimated phenology metrics compared are as follows: start-of-season time (SOST), start-of-season NDVI (SOSN), end-of-season time (EOST), end-of-season NDVI (EOSN), maximum time (MAXT), and maximum NDVI (MAXN). The comparison statistical results are the five-year average (2017–2021) coefficient of determination (R²), Pearson’s r (r), mean absolute error (MAE), and relative root mean square error (RRMSE). The number of samples for cheatgrass was 7813, and medusahead had 3693 total points.

		SOST		SOSN		EOST		EOSN		MAXT		MAXN
		BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8
2017–2021	R2	0.95	0.97	0.91	0.91	0.84	0.53	0.93	0.86	0.87	0.88	1.00	1.00
	r	0.98	0.98	0.95	0.95	0.92	0.73	0.97	0.93	0.94	0.94	1.00	1.00
	MAE	0.19	0.11	0.22	0.11	0.44	0.29	0.69	0.56	0.37	0.23	0.00	0.00
	RRMSE	4.73	2.82	0.56	0.35	3.75	3.43	1.40	1.36	3.26	2.32	0.02	0.05

). Given fewer samples for modeling compared to cheatgrass, the medusahead regression tree did well at estimating phenological metrics relative to the phenology decision tree results. EOST had the lowest correlation with an average Pearson’s r of 0.73 (R² = 0.53). Extracting the results from the medusahead estimations and comparing them to the 20% withheld from the model (Table 2) revealed that the maps had a high correlation with the temporal metrics and maintained a strong relation with NDVI metrics.

Table A1. Cheatgrass (BRTE [44]) and medusahead (TACA8 [44]) estimated phenology metrics compared to decision-tree-produced phenology. The estimated phenology metrics compared are as follows: start-of-season time (SOST), start-of-season NDVI (SOSN), end-of-season time (EOST), end-of-season NDVI (EOSN), maximum time (MAXT), and maximum NDVI (MAXN). The comparison statistical results are the coefficient of determination (R²), Pearson’s r (r), mean absolute error (MAE), and relative root mean square error (RRMSE).

		SOST		SOSN		EOST		EOSN		MAXT		MAXN
		BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8	BRTE	TACA8
2017	R2	0.96	0.98	0.94	0.92	0.71	0.55	0.86	0.84	0.78	0.90	1.00	1.00
	r	0.98	0.99	0.97	0.96	0.86	0.75	0.94	0.92	0.88	0.95	1.00	1.00
	MAE	0.21	0.11	0.24	0.12	0.76	0.57	1.35	0.86	0.45	0.34	0.00	0.00
	RRMSE	5.55	2.76	0.57	0.36	5.53	5.22	2.27	1.79	3.72	2.90	0.01	0.05
2018	R2	0.96	0.97	0.91	0.95	0.87	0.61	0.94	0.96	0.90	0.85	1.00	1.00
	r	0.98	0.98	0.95	0.98	0.94	0.78	0.97	0.98	0.95	0.92	1.00	1.00
	MAE	0.13	0.05	0.17	0.03	0.42	0.12	0.64	0.31	0.32	0.16	0.00	0.02
	RRMSE	3.84	2.25	0.60	0.22	3.58	2.28	1.36	0.97	2.97	1.93	0.02	0.17
2019	R2	0.98	0.99	0.95	0.98	0.86	0.41	0.94	0.91	0.92	0.89	1.00	1.00
	r	0.99	1.00	0.98	0.99	0.94	0.64	0.97	0.96	0.96	0.95	1.00	1.00
	MAE	0.11	0.04	0.13	0.04	0.40	0.19	0.54	0.36	0.31	0.15	0.00	0.00
	RRMSE	2.90	1.31	0.38	0.18	3.57	3.10	1.20	1.14	2.88	1.88	0.02	0.03
2020	R2	0.96	0.99	0.94	0.95	0.90	0.55	0.96	0.90	0.91	0.92	1.00	1.00
	r	0.98	0.99	0.97	0.97	0.95	0.75	0.98	0.95	0.95	0.96	1.00	1.00
	MAE	0.16	0.08	0.16	0.09	0.32	0.26	0.45	0.46	0.32	0.18	0.00	0.00
	RRMSE	4.53	2.00	0.42	0.28	2.98	3.41	1.03	1.24	2.97	2.00	0.00	0.00
2021	R2	0.91	0.92	0.80	0.75	0.83	0.56	0.95	0.72	0.87	0.86	1.00	0.98
	r	0.96	0.96	0.90	0.87	0.91	0.75	0.98	0.87	0.94	0.94	1.00	0.99
	MAE	0.33	0.29	0.38	0.29	0.29	0.28	0.48	0.80	0.44	0.32	0.00	0.00
	RRMSE	6.81	5.76	0.82	0.70	3.09	3.15	1.13	1.65	3.77	2.87	0.04	0.00

displays each year’s estimated phenology compared to the decision tree analysis.

The estimated phenology maps were trained on cheatgrass- and medusahead-dominated areas but mapped across all vegetation types in the SRP and NBR. We evaluated 48 predictor variables for the cheatgrass models and the medusahead models, and the top variables were noted as most influential across the model tree in Figure 7. The SOST for cheatgrass was driven greatly by week 14 NDVI values and medusahead by week 15 NDVI. Changes in those top influential drivers can impact the SOST values by an increase or decrease in NDVI values. The remainder of the variables had lower relative driver factors but collectively influenced the outcome. EOST showed more variance in usage throughout all predictor variables, which is why the “Other” category has a larger percentage than its counterparts in Figure 6. This showed the difficulty in identifying EOST with active vegetation present and not being able to rely on a few NDVI variables.

The phenology metric maps displayed the differences between the start of sustained active growth of cheatgrass and medusahead throughout the entire study area. Figure 8 illustrates the final estimated maps throughout the six years mapped. Pixels labeled with at least 1% medusahead more than doubled in 2022 compared to previous years, with a majority of the increase between 1% and 5% cover [20]. An earlier SOST for cheatgrass (SOST less than week 12) and medusahead (SOST less than week 14) was noticed in lower-elevation areas in the SRP ecoregion. These earlier SOST areas also spatially aligned with areas of at least 15% cheatgrass or 5% medusahead cover in the six-year (2017–2022) average. Figure A1 illustrates the annual cover for cheatgrass (Figure A1a–f) and medusahead (Figure A1g–l). Phenology results for 2022 SOST (Figure 8f and Figure 8l, respectively) and MAXT captured similar growing season lag for medusahead’s phenology compared to cheatgrass’s in the previous years. Data for 2022 maps were not used to train phenology models, so identifying similar latency patterns in SOST and MAXT with pixels of at least 15% cover shows some robustness in these models. These models continued to show robustness by differentiating cheatgrass and medusahead phenology based on the same NDVI time series unseen by the models.

Using the stratified pixels described in Section 2.3 and Figure 3a, we extracted the values from the phenology maps for comparison (104,963 points for cheatgrass and 52,070 points for medusahead). Figure 9 shows the estimated values of medusahead SOST for the years 2017–2022, which had an average difference of +2.7 weeks compared to cheatgrass. The MAXT showed a closer difference between medusahead and cheatgrass with an average of 2 weeks from 2017 to 2022. The Figure 9 distribution showed similar ranges to those shown by training data in Figure 4. The estimated data values’ interquartile ranges increased in Figure 9, capturing more variety within the larger sample size for medusahead’s MAXT, displaying the delayed greenness as compared to cheatgrass. Figure 9 captures the annual SOST and MAXT differences for cheatgrass and medusahead and shows a reduction in differences throughout the years. Based on the NDVI profiles for pixels with 20% cover, the annual NDVI values are also decreasing with the species. We found that the extracted average week 14 and week 15 NDVI also had these similar differences between cheatgrass and medusahead. In Figure 9, the reduction in the difference between medusahead and cheatgrass was noticed for SOST and MAXT throughout the years because of the decrease in NDVI values for these models’ most important drivers (week 14 and week 15 NDVI (Figure 7)). Training data values for 2021 (Figure 3) showed similar temporal patterns and ranges for cheatgrass and medusahead, but medusahead’s delayed greenness compared to cheatgrass was still captured in the maps (Figure 8). Capturing these different growing patterns was important because this was a notable difference between cheatgrass and medusahead phenological patterns. The models were able to annually capture the delayed latent growing seasons of medusahead compared to cheatgrass while using the same NDVI datasets for both species’ models.

Comparing cheatgrass and medusahead phenology metrics’ median values, the combined years showed all temporal metrics for medusahead had a significantly later seasonal pattern than cheatgrass (

Table A2. Cheatgrass and medusahead Mann–Whitney U-test statistics for 2017–2022. U-statistic with a p-value < 0.01 is indicated by “**”, 0.01 < p-value < 0.05 is indicated by “*”, and p-value > 0.05 is indicated by “ ”. Effect size was calculated by rank-biserial correlation; a value of 0 indicates that cheatgrass and medusahead were present in equal proportions, a value of −1 indicates complete coverage by medusahead, and a value of 1 indicates complete coverage by cheatgrass. MAXN 2022 p-value = 0.28.

		SOST	SOSN	EOST	EOSN	MAXT	MAXN
2017	U-statistic	197,468,344.5 **	89,387,548.5 **	141,249,540.5 **	85,685,862.5 **	187,547,484.5 **	35,507,098 **
	effect size	−0.85	0.16	−0.32	0.20	−0.76	0.67
2018	U-statistic	210,365,165 **	93,046,952.5 **	151,436,039 **	113,544,012.5 **	197,023,902.5 **	75,791,524 **
	effect size	−0.82	0.20	−0.31	0.02	−0.70	0.35
2019	U-statistic	180,185,628.5 **	146,370,580 **	188,988,769 **	182,153,827 **	215,571,194.5 **	164,276,548 **
	effect size	−0.40	−0.14	−0.47	−0.42	−0.68	−0.28
2020	U-statistic	147,928,271.5 **	125,572,657 **	179,987,034.5 **	176,333,033 **	190,429,299.5 **	163,344,940.5 **
	effect size	−0.38	−0.17	−0.68	−0.65	−0.78	−0.53
2021	U-statistic	56,210,773 **	68,597,582 **	70,585,046 **	74,527,591 **	59,946,545 **	77,348,693.5 **
	effect size	−0.33	−0.62	−0.66	−0.76	−0.41	−0.82
2022 a	U-statistic	10,676,011.5 **	10,289,900 *	14,793,177 **	7,371,812 **	13,857,482.5 **	9,785,941
	effect size	−0.07	−0.04	−0.49	0.26	−0.39	0.02
2017–2022	U-statistic	4,475,995,305 **	3,157,955,127.5 **	4,320,763,126.5 **	3,824,070,167.5 **	4,862,155,702.5 **	3,256,081,399 **
	effect size	−0.56	−0.10	−0.51	−0.34	−0.70	−0.14

^a SOSN 2022 p-value = 0.012 and MAXN 2022 p-value = 0.28.

). Some of the individual years were not as strongly significant with a rank-biserial correlation of closer to 0 than −0.4. The overall largest significant difference was found for MAXT, where the combined years had an effect size of 70% of medusahead values being larger than cheatgrass.

Figure 10 violin plots show the variance distribution of the NDVI values for SOSN and MAXN around the median values. The wider the area on the plot, the more densely populated that NDVI value. NDVI metrics had a different pattern emerge because the NDVI metrics were extracted based on the modeled temporal values. Although MAXT was significantly different (p < 0.05), MAXN did not show this same significant difference. The combined years for MAXN had an effect size of only 14% of medusahead values being greater than cheatgrass (Table A2). The only comparison that had no significant difference (p > 0.05) was MAXN for 2022 (p = 0.28). This was likely due to the pixels evaluated having similar cover estimates.

4. Discussion

We were able to differentiate medusahead and cheatgrass phenology patterns at 30 m resolution by training a regression tree phenology model on pixels with at least 20% cover of cheatgrass or medusahead. Because the SOST, EOST, and MAXT test data were significantly (p < 0.05) different between the species, we infer that based on our phenology metrics we can differentiate phenological differences between cheatgrass and medusahead at 30 m spatial resolution. The temporal metrics were highly significant with a large effect size in favor of medusahead having a later growing season than cheatgrass. The overall insignificance of MAXN for the combined years was expected because the relative maximum between the SOST and EOST was close to the same NDVI value for both species. Because cheatgrass and medusahead have similar vegetation behavior and structure, NDVI values are likely similar. Our results indicate that SOSN, EOSN, and MAXN values can be used to identify the NDVI thresholds for cheatgrass and medusahead in the SRP and NBR ecoregions.

To interpret the NDVI values, we focused on the sustained growth of cheatgrass and medusahead. The sustained growth was interpreted as the 5 weeks of sequential increasing scNDVI. Although the NDVI pixels have a heterogeneous community of vegetation, the increase in one scNDVI value is equivalent to one one-hundredth of an increase in raw NDVI, which is a relatively large increase compared to raw NDVI. The differences in sustained growth between species we found agree with previous studies that observed two-to-three-week divergences between vegetation stages (i.e., heading and flowering stages) dates for these two species in Idaho [13,52]. Generally, NDVI values increase as vegetation progresses through life cycles up to the flowering stage, in which the maximum NDVI is achieved, and then decreases as the plants senesce [53]. For cheatgrass and medusahead, the flowering stage occurs after the boot stage, and this transitional period corresponds to MAXT, which aligned with observations conducted in the Idaho area [13,52]. The flowering stages for cheatgrass and medusahead were reported as weeks 21 and 24 [13], respectively, and 21 and 23, respectively [52]. In our study, for years 2017–2022, the median estimated cheatgrass and medusahead MAXT ranged from week 19 to 21 (cheatgrass) for years 2017–2022, and week 20 to 23 (medusahead). This timing aligned with observed flowering stages. The MAXT maps showed similar results to those findings from year to year.

Several studies have shown that the growing seasons for cheatgrass and medusahead are different [13,18,24,52], and that difference was observed in our training data. The training data had median phenology values that averaged 3.5 weeks difference for SOST, and 2.4 weeks difference for MAXT between cheatgrass and medusahead. These differences boosted the authors’ confidence in using these training datasets for developing phenology models. Developing the decision tree rules involved several steps for reiterating the manual processing for identifying active growth. Developing these decision tree rules reduced human error and generated optimized phenometrics, which lowered errors. This approach may, however, cause some model errors that human eyes could detect, such as identifying if an area has suitable EAG cover or an area that was not masked out by NLCD [21]. To prevent these errors in the training data, we were able to modify restrictions and omit potential outliers. We then produced model-estimated annual phenology metrics that were highly correlated with the manually influenced automated decision tree phenology that was withheld from the training data.

Extracting MAXT for NDVI is relatively straightforward within an NDVI time series of a monoculture. However, given exotic annual grasses often grow intermixed with perennial forbs and other rangeland species that have a similar growing season [18,22], extracting maximum time for NDVI is difficult. A heterogeneous vegetation community within 30 m pixels makes extracting the maximum NDVI values for a particular species challenging [38]. Previously, we extracted the local maximum NDVI values which created uniform maximum metrics between cheatgrass and medusahead [51]. We opted for modeling cheatgrass and medusahead MAXT, which extracted the delayed maximum growth between cheatgrass and medusahead based on the same NDVI time series given the desired pixel. Pixels with high medusahead and cheatgrass cover (≥20%) produced similar temporal values (SOST comparison r² = 0.97 and MAXT comparison r² = 0.86). The mixed pixel effect decreased the temporal estimates compared to their respective medians in Figure 9 for medusahead in the years 2017–2022 by one week and increased cheatgrass by one week. The mixed pixel effect influenced the phenology temporal metrics (by ±1 week) between pixels with the same amount of species cover (at least 20% cover) [38].

One pattern emerged in the phenology metric maps where pixels with higher cheatgrass or medusahead cover had an earlier start of sustained growth. This pattern could be an indication of capturing earlier sustained growth by the dominant species in the pixel [38]. When in the field, we observed a plot with early season medusahead emergence covered by the previous year’s medusahead plant residue. This plot was south of our study area and showed useful growing observations of medusahead (e.g., emerging within dormant species or blanketing the ground with last year’s growth) which gave the impression of early emergence near Reno, Nevada (during a visit on 31 January 2024). At the end of the week, snow fell and covered the plot after a weekly average maximum temperature of 15.2 °C [54], which could have halted the detectable growth of medusahead or suppressed its reflectance. This early emergence seems plausible because the average temperature for the visiting week was above the average normal maximum temperature of 8.6 °C (measured from 1981 to 2010 for 28 January–2 February) [55]. An earlier onset of spring thaw can occur with warmer than average temperatures and has been indicated as a driver for cheatgrass and medusahead growth [8,18,19]. The pixels with more cheatgrass or medusahead cover should have a more homogeneous reflectance cover because the photosynthetic activity will be similar, making it easier to detect sustained growth, and indicate earlier greenness with less influence of dormant native vegetation’s NDVI values [38]. The NDVI values of pixels with lower cover of the species could be influenced by the dormant vegetation that stands above the emerging greenness, suppressing the plant reflectance early in the season. This may affect the uniformity throughout the maps where the same cover values differ in temporal signatures, but more data filtering would be needed to prevent mixed pixel effects that may influence phenology patterns [38].

We focused on the SOST and MAXT during the sustained growing period because these metrics have been shown to influence timing for mitigation [25,26]. Herbicide applications have reduced cheatgrass and medusahead cover during vegetation and reproductive stages and prior to maturity. Herbicide applications can be planned using SOST and MAXT [26]. SOST and MAXT can also be used to plan grazing. Research has found greater nutritional values prior to MAXT for medusahead [25]. Mowing was also shown to be more effective at reducing medusahead cover before the final stages of seed development (past the peak portion of the growing season) according to Brownsey et al. [25]. Our annual phenology products for cheatgrass and medusahead can be used to inform the planning of those mitigation approaches.

5. Conclusions

Overall, we found developing species-oriented models for cheatgrass and medusahead to be challenging. Although difficult to produce, species-oriented models were used to identify past active growing seasons of sustained growth. Developed maps can be used to find patterns among the differences in these species. Developing training data for additional years would create a more robust model that could be applied to a broader study area, potentially expanding the utility of the modeling effort.

Developing each metric model individually using its own set of parameters and hyperparameters could help improve future phenology estimates by further constraining the modeled variables to features within the metrics’ range. For example, when SOST did not surpass week 23 or MAXT only occurred between a certain range of weeks, the model could omit NDVI variables later than week 23 or not within the specified range. These restrictions could help focus the model on the more key features for development. Another consideration for increasing accuracy is to loosen EOST restrictions based on specific species, such as the continued growth of medusahead after cheatgrass matures.

These phenology metrics can be used by rangeland managers to inform planning treatments designed to reduce the spread of these exotic species. The model-estimated results include SOST, MAXT, and EOST, followed by pixel drilling to extract SOSN, MAXN, and EOSN, and calculated DUR and AMP. Aligning herbicide or grazing timing with the timing of SOST and MAXT may best target the annual grasses at peak vulnerability. Phenology characteristics can be used to inform the development of efficient management plans by identifying optimal peak timing to control cheatgrass and medusahead with less effort.