Changes in Meadow Phenology in Response to Grazing Management at Multiple Scales of Measurement

: Riparian and ground-water dependent ecosystems found in the Great Basin of North America are heavily utilized by livestock and wildlife throughout the year. Due to this constant pressure, grazing can be a major inﬂuence on many groundwater dependent resources. It is important for land managers to understand how intensity and timing of grazing affect the temporal availability of these commodities (i.e., biodiversity, water ﬁltration, forage, habitat). Shifts in forage or water availability could potentially be harmful for fauna that rely on them at speciﬁc times of the year. Seven meadow communities, each consisting of three distinct vegetative communities, were grazed at three intensities to determine the relationship between grazing management and phenological timing of vegetation. The agreement of on-the-ground measurements, near-surface digital cameras (phenocams), and satellite-based indices of greenness was examined for a two-year period (2019–2020) over these grazing and vegetative community gradients. Field determined phenology, phenocam Green Chromatic Coordinate ( GCC ), and Landsat Normalized Difference Vegetation Index ( NDVI ) were all highly correlated and the relationship did not change across the treatments. Timing of growth varied in these ecosystems depending on yearly precipitation and vegetative type. Communities dominated by mesic sedges had growing seasons which stopped earlier in the year. Heavier grazing regimes, however, did not equate to signiﬁcant changes in growing season. Ultimately, shifts in phenology occurred and were successfully monitored at various spatial and temporal scales. and uncontrolled meadows. It should be noted that in 2019, all meadows were treated as uncontrolled. In 2020, grazing treatment actually began. Designating meadows as grazed, ungrazed, and uncontrolled in 2019 was done solely to show differences in the groups of meadows across the years.


Introduction
Understanding how groundwater dependent areas located within rangelands are affected by various management techniques is crucial to conservation. Their importance in arid and semi-arid ecosystems is generally disproportionate to their overall size [1,2]. Riparian areas boost landscape biodiversity [3], filter water [4], provide forage and water for livestock [5], and serve as suitable seasonal and year-round habitat for a wide variety of key taxa [6,7].
Livestock tend to congregate in these systems, due to increased water and forage availability [8]. Improper grazing management has been shown to reduce vegetation, leading to destabilized streambanks and erosion [9,10]. Over-grazing can also decrease leaf area, nutrient stocks, photosynthetic capability, reproductive success, and vegetative growth [11][12][13]. It should be noted that varying plant communities can lead to different responses in grazing. Skaer et al. [14] found that grazing at Palo Corona Regional Park (Monterey County, California) led to a 15% increase in cover of exotic annual forbs but had little to no effect on the cover of native perennial forbs and grasses.
One of the main obstacles that land managers face is understanding how disturbance affects the high temporal variability that accompanies plant establishment and growth [15].
GCC, and Landsat NDVI data within these heterogeneous systems and across a spectrum of different grazing regimes.

Site Description and Treatments
Several high elevation meadows, referred to as the Haypress meadow complex, are located in the Desatoya Mountains of the Central Great Basin (39 • 27 N, 117 • 36 W; 39 • 19 N, 117 • 42 W). Over the past 10 years (2010-2020), annual precipitation in this complex ranged between 175 and 419 mm, with an average value of 315 mm. Annual mean temperature ranged from 4.6 to 7.6 • C, with an average of 5.8 • C. These values were obtained via the Google Earth Climate Engine (http://climateengine.org (accessed on 1 October 2021), [44]). Four meadows were selected for this study; they were given names based on their relative locations or the grazing intensity associated with them: (1) Lower at 2322 m, (2) Middle at 2365 m, (3) Upper at 2438 m, and (4) Uncontrolled at 2320 m. For the Lower, Middle, and Uncontrolled meadows the soil type for the first 14 inches of soil was a gravelly loam; lower than 14 inches the soil type was a gravelly clay loam. In the Upper meadow, the soil type for the first 12 inches was a gravelly clay loam; beneath that it became an extremely cobbly loam [45].
The vegetation of the meadows, while diverse, was delineated into three distinct classes. These classes were (1) the dry meadow community, which consisted of an overstory of mountain big sagebrush (Artemisia tridentata ssp. vaseyana) and an understory of primarily Douglas' sedge (Carex douglasii), (2) the mesic meadow community, which consisted of Douglas' sedge, arctic rush (Juncus arcticus), and a mixture of meadow forbs (ex. Chorispora tenella, Symphyotrichum ascendens), and (3) the wet meadow community, which consisted of Nebraska sedge (Carex nebrascensis), arctic rush, and a mixture of facultative wet forbs (ex. Montia chamissoi, Veronica americana). Vegetation surveys were conducted in 2020, using standard line-point intercept methods [46], to verify the differentiation of these groups (Figure 1). On average, forbs had the highest amount of cover in the wet communities (29%), graminoids in the mesic communities (71%), and shrubs in the dry communities (56%).
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 21 upland meadows within the Great Basin, (2) quantify the effects of differing grazing management strategies on phenological metrics in the context of year-to-year climate variability, and (3) determine the relationship between on-the-ground phenology measurements, phenocam GCC, and Landsat NDVI data within these heterogeneous systems and across a spectrum of different grazing regimes.

Site Description and Treatments
Several high elevation meadows, referred to as the Haypress meadow complex, are located in the Desatoya Mountains of the Central Great Basin (39°27' N, 117°36' W; 39°19' N, 117°42' W). Over the past 10 years (2010-2020), annual precipitation in this complex ranged between 175 and 419 mm, with an average value of 315 mm. Annual mean temperature ranged from 4.6 to 7.6°C, with an average of 5.8°C. These values were obtained via the Google Earth Climate Engine (http://climateengine.org, [44]). Four meadows were selected for this study; they were given names based on their relative locations or the grazing intensity associated with them: (1) Lower at 2322 m, (2) Middle at 2365 m, (3) Upper at 2438 m, and (4) Uncontrolled at 2320 m. For the Lower, Middle, and Uncontrolled meadows the soil type for the first 14 inches of soil was a gravelly loam; lower than 14 inches the soil type was a gravelly clay loam. In the Upper meadow, the soil type for the first 12 inches was a gravelly clay loam; beneath that it became an extremely cobbly loam [45].
The vegetation of the meadows, while diverse, was delineated into three distinct classes. These classes were 1) the dry meadow community, which consisted of an overstory of mountain big sagebrush (Artemisia tridentata ssp. vaseyana) and an understory of primarily Douglas' sedge (Carex douglasii), 2) the mesic meadow community, which consisted of Douglas' sedge, arctic rush (Juncus arcticus), and a mixture of meadow forbs (ex. Chorispora tenella, Symphyotrichum ascendens), and 3) the wet meadow community, which consisted of Nebraska sedge (Carex nebrascensis), arctic rush, and a mixture of facultative wet forbs (ex. Montia chamissoi, Veronica americana). Vegetation surveys were conducted in 2020, using standard line-point intercept methods [46], to verify the differentiation of these groups ( Figure 1). On average, forbs had the highest amount of cover in the wet communities (29%), graminoids in the mesic communities (71%), and shrubs in the dry communities (56%). Historically, grazing in these meadows consisted of yearlong grazing by feral horses and managed cattle grazing. Haypress livestock grazing occurred in the following manner: for two consecutive years, grazing occurred in May and continued through mid-June; the following two years, use occurred mid-August through September. Livestock use was generally limited to one month, with around 700 cow-calf pairs having access to the meadows but could last up to six weeks depending on the year. In 2019, cattle were brought in on the first half of the grazing rotation (mid-June).
In the winter of 2019, pipe-rail fencing was installed in the complex. The perimeters of the Lower, Middle, and Upper meadows were fenced completely, while the Uncontrolled meadow was not fenced. For the three fenced meadows, cross fencing was constructed which cut each meadow in half. The higher elevation portion of each meadow was designated as part A, and the lower part B. This allowed for a total of seven separate units to be treated and measured (Table 1). In the summer of 2020, cattle were brought into the complex in mid-June. Section A of each meadow was grazed such that the major forage species of the riparian areas had a mean stubble height of roughly 10 cm [47] (approximately 3 days of grazing). Section B of each meadow was left ungrazed, and the Uncontrolled meadow was left unmanaged (allowed for constant grazing by both cattle and horses to occur in the meadow throughout the summer).

Field Methodology
The Line-Point Intercept (LPI) method [46] was used to calculate the percent cover of the wet, mesic, and dry subregions of each of the seven meadow sections in June of 2018. The number of transects, length of transects, and number of points along transects was tailored to each area of vegetation ( Figure 2, lengths of transects ranged from 3 to 20 m, number of transects per community ranged from 3 to 9). Monitoring was set up in such a way that at least 100 data points could be generated from each LPI measurement (points ranged from 100 to 120).
In 2019 and 2020, these same transects were used to monitor the phenology and height of the meadows bimonthly, throughout the growing season (May to October). At 25%, 50%, and 75% of the total length of each transect, a 1 m hoop was placed and average height (collected using a fiberglass folding ruler) and phenology of each species in the hoop was recorded ( Figure 2). Average phenology was determined using simplified methods outlined in Moore et al. [48], where each species was assigned one of 8 different phenological stages: Leaf development (L), Initial inflorescence (I), Mature inflorescence (Im), Inflorescencepollen released (Ip), Initial seed development (S), Mature seed development (Sm), Dropped seed (Sd), and Desiccated (D). If a species was grazed to the point where phenology could not be ascertained, a rating of G (Grazed) was given. A total of roughly 9-15 hoops were placed per community type for each of the meadow regions. In 2019 and 2020, these same transects were used to monitor the phenology and height of the meadows bimonthly, throughout the growing season (May to October). At 25%, 50%, and 75% of the total length of each transect, a 1 m hoop was placed and average height (collected using a fiberglass folding ruler) and phenology of each species in the hoop was recorded ( Figure 2). Average phenology was determined using simplified methods outlined in Moore et al. [48], where each species was assigned one of 8 different phenological stages: Leaf development (L), Initial inflorescence (I), Mature inflorescence (Im), Inflorescence-pollen released (Ip), Initial seed development (S), Mature seed development (Sm), Dropped seed (Sd), and Desiccated (D). If a species was grazed to the point where phenology could not be ascertained, a rating of G (Grazed) was given. A total of roughly 9-15 hoops were placed per community type for each of the meadow regions.

Phenocam and Landsat Methodology
In fall of 2017, four phenocams (StarDot NetCam SC 5MP IR-enabled cameras using complementary metal oxide semiconductor (CMOS) image sensors, StarDot Technologies, Buena Park, CA, USA) were installed in four of the meadow sections previously discussed, in fall of 2018 the other three were installed. Methods of installation and use for the phenocams were sourced from the PhenoCam Network

Phenocam and Landsat Methodology
In fall of 2017, four phenocams (StarDot NetCam SC 5MP IR-enabled cameras using complementary metal oxide semiconductor (CMOS) image sensors, StarDot Technologies, Buena Park, CA, USA) were installed in four of the meadow sections previously discussed, in fall of 2018 the other three were installed. Methods of installation and use for the phenocams were sourced from the PhenoCam Network (https: //phenocam.sr.unh.edu/webcam/tools/ (accessed on 1 October 2021)) ( Table 1). Additionally, time-domain reflectometry soil moisture sensors (CS-616, Campbell Scientific, Logan, UT, USA) were installed at depths of 10, 20, 35, and 50 cm at each of the phenocam sites. Eight RGB images at a resolution of 1296 × 960 pixels were taken daily between 12:00 pm and 3:30 pm. All data and images were downloaded monthly during the growing season, and once every three months during the winter.
Phenopix, the R statistical program package [49], was used for the calculation of all phenological metrics derived from the phenocam images. The package allowed Regions of Interest (ROIs) to be manually defined that were representative of the wet, mesic, and dry communities within each meadow section. The average digital numbers, which ranged from 0 to 255 and were extracted for red (RDN), green (GDN), and blue (BDN), of the pixels Remote Sens. 2021, 13, 4028 6 of 19 within each ROI were used for the analysis. These numbers were converted into green chromatic coordinate (GCC) values, a relative percent index, using the following equation: GCC values were filtered using the night, spline, and max filter options in the phenopix package. An additional novel filter using the blue chromatic coordinate, was added to remove values influenced by snow (see Snyder et al. [43] for more details). A method proposed by Gu et al. [50] was used to apply a double logistic fit to the filtered sub daily GCC values. One thousand replications were used for an uncertainty analysis to determine how well predicted values fit observed. The Gu et al. [50] threshold method calculates four phenophase dates for each fitted set of data: an upturn date (UD) when GCC of vegetation begins to increase consistently, a stabilization date (SD) when vegetation approaches maximum GCC, a downturn date (DD) when GCC starts to consistently diminish, and a recession date (RD) when vegetation reaches a seasonal low. The double logistic fit determined by Gu et al. [46] was also applied to the collected soil moisture data and the RD phenophase DOY was determined for each individual meadow across the two years.
Landsat NDVI data were extracted from Landsat 8 Optical Imager (OLI) satellite images using the Climate Engine (http://climateengine.org (accessed on 1 October 2021)) [44]. Polygons were created for each meadow section and NDVI was determined for 2019 and 2020. NDVI values were calculated as Equation (2): where NIR is the near infrared at-surface reflectance and Red is the red at-surface reflectance. The polygons created were large enough to encompass dry, mesic, and wet vegetation communities, and were at least one pixel in size (30 m × 30 m). Climate Engine applied a cloud mask to the Landsat top of atmosphere (TOA) reflectance data. The cloud masking removed medium and high confidence snow, shadow, and cirrus clouds using the Quality Assessment Band provided in the Landsat Google Earth Engine (GEE) collection.

Statistical Methodology
A random forest model was used to determine which variables had the greatest importance in predicting the on-the-ground phenological state of vegetation. Random forest models are used in the ecological community when dealing with issues such as auto-correlation [51,52]. This is due to these models being non-parametric, using bootstrap aggregation, and implementing feature randomness. Table 2 shows the variables used in the analysis. The field data collected was input to the randomForest package [53] in the R statistical program [54] to create and test the model, with four simplified phenological stages used as the output: Leaf (L), Inflorescence (I), Seed (S), and Desiccated (D). The analysis was completed using 1000 trees and was set to test three variables at each split. Out-of-box estimate of error was used to measure the strength of the model and the mean decrease in Gini coefficient was used to measure the importance of each variable. To further test the strength of correlation between field-and camera-based observations, the relationship between the four field-based, simplified phenological stages and the four camera-based phenophase dates was analyzed. The date where the highest percentage of each simplified phenological stage was seen in each community type of each meadow region was calculated. These values were then plotted against their corresponding camerabased phenophase dates (UD-L, SD-I, DD-S, and RD-D) and repeated measures correlation coefficients were calculated, with meadows used as the subject for repeated measure. Factorial ANOVAs and Tukey tests were utilized to determine if there were differences in UD, SD, DD, and RD dates across vegetation communities, years, and grazing interactions. Meadow location was used as a randomized block variable. Since the grazing treatment began in 2020, the analysis was treated as a before-after design where a significant grazing treatment effect would appear in the grazing-year interaction. Similar ANOVAs and Tukey tests were created for field phenology dates (max observance of L, I, S, and D). Data met the assumptions of homoscedasticity, independence, and normality. Root mean square error (RMSE) and Pearson's r were used to determine the level of correlation between Landsat derived NDVI and camera-derived GCC. Comparisons between Landsat NDVI and camera-derived GCC were completed using data from the dates of the Landsat imagery.

Phenology Stage Modeling
Grazing treatments were successfully implemented and monitored in 2019 and 2020. In 2020, grazed meadows were brought down to a 10 cm stubble height in the wet communities, with cattle tending to focus their attention on the wetter areas and leaving the other sections of the meadows largely undisturbed for the period they were allowed in the exclosures. Phenological data of the vegetation across all meadows were collected, both in person and remotely, and a random forest model was successfully created to predict phenological stages in the Haypress complex (L-Leaf; I-Inflorescence; S-Seed; D-Desiccated). The out-of-box estimate of error (OOBE) obtained by the model was 33.56%. Across the phenological classes, there were significant differences in error. Phenological stages D and L were correctly predicted most frequently by the model, with class errors of the stages being 0.2772 and 0.1985, respectively. The model had a more difficult time with predicting stages I and S, with class errors of 0.5962 and 0.8703. DOY, GCC, life form, and soil moisture were the variables with the highest level of importance based on their mean decreases in Gini. Camera, year, grazing, and meadow were the variables with the lowest (Figure 3).
Phenological stages D and L were correctly predicted most frequently by the model, with class errors of the stages being 0.2772 and 0.1985, respectively. The model had a more difficult time with predicting stages I and S, with class errors of 0.5962 and 0.8703. DOY, GCC, life form, and soil moisture were the variables with the highest level of importance based on their mean decreases in Gini. Camera, year, grazing, and meadow were the variables with the lowest (Figure 3).

Phenocam and Field-Based Metrics and Correlations
Different vegetation communities were effectively delineated from the phenocam images and phenophase dates were extracted (UD-upturn date; SD-stabilization date; DD-downturn date; RD-recession date) (Figure 4).

Phenocam and Field-Based Metrics and Correlations
Different vegetation communities were effectively delineated from the phenocam images and phenophase dates were extracted (UD-upturn date; SD-stabilization date; DD-downturn date; RD-recession date) (Figure 4).
A statistically significant difference in DOY for SD, DD, and RD phenophases was found across both vegetation community type and year, however not for meadow or grazing-year interaction. For the UD phenophase, significant differences were only seen across years ( Table 3). The mean difference between 2019 and 2020 was highest for the DD phenophase (−42.13) and lowest for the UD phenophase (−14.13). A Tukey post-hoc test revealed significant pairwise differences only for the SD mesic-dry comparison in 2019. In 2020, community differences were seen for mesic-dry and mesic-wet comparisons across all phenophases except UD ( Figure 5). No differences between wet-dry communities were statistically significant in either 2019 or 2020.
The date when each phenological stage was most frequently observed in the field was unique for both years. Consistently, across both meadows and years, the phenology stages that had the closest maximum observed DOYs were I and S ( Figure 6). Significance of variables for field-observed data was similar to that previously discussed for phenocam data. Meadow and grazing-year interaction were not significant for any phenology stage, while year was significant for all stages (Table 4). Vegetative community was only significant for S, with significant interannual differences solely occurring between the mesic and other two communities in 2020. Correlation between camera-based phenophase dates and field-based phenological stages was very high (Figure 7). Year-to-year changes led to some variability in the repeated measures correlation coefficients of the data, with a range of 0.725-0.942 across the two years and all community types. All correlation values were above 0.90, except for the mesic in 2020 (0.725). Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 21 A statistically significant difference in DOY for SD, DD, and RD phenophases was found across both vegetation community type and year, however not for meadow or grazing-year interaction. For the UD phenophase, significant differences were only seen across years ( Table 3). The mean difference between 2019 and 2020 was highest for the DD phenophase (−42.13) and lowest for the UD phenophase (−14.13). A Tukey post-hoc test revealed significant pairwise differences only for the SD mesic-dry comparison in 2019. In 2020, community differences were seen for mesic-dry and mesic-wet comparisons across all phenophases except UD ( Figure 5). No differences between wet-dry communities were statistically significant in either 2019 or 2020.    The date when each phenological stage was most frequently observed in the field was unique for both years. Consistently, across both meadows and years, the phenology     Grazing intensity, as shown in the ANOVA model, had little effect on phenophase date (Figure 8). For the UD, SD, and DD phases, the treatment that had the latest average date varied across community type. For example, for the UD phase in 2020, the grazed Grazing intensity, as shown in the ANOVA model, had little effect on phenophase date (Figure 8). For the UD, SD, and DD phases, the treatment that had the latest average date varied across community type. For example, for the UD phase in 2020, the grazed meadows had the latest average date in the dry communities (102), ungrazed had the latest average date in the mesic communities (109), and the uncontrolled meadow had the latest date in the wet communities when compared to the averages of the other treatment groups (110). Additionally, when looking at the same community type, the treatment with the latest average date was never the same across all phenophases. The RD phase was slightly different in that, in 2019, the meadow that would be uncontrolled had the latest date in the dry communities, and the meadows which would be ungrazed had the latest average date for the wet and mesic communities. In 2020, however, the uncontrolled meadow had the latest date compared to the averages of the other treatments across all communities. It should be noted that these shifts did not equate to significance in the ANOVA model.  (109), and the uncontrolled meadow had the latest date in the wet communities when compared to the averages of the other treatment groups (110). Additionally, when looking at the same community type, the treatment with the latest average date was never the same across all phenophases. The RD phase was slightly different in that, in 2019, the meadow that would be uncontrolled had the latest date in the dry communities, and the meadows which would be ungrazed had the latest average date for the wet and mesic communities. In 2020, however, the uncontrolled meadow had the latest date compared to the averages of the other treatments across all communities. It should be noted that these shifts did not equate to significance in the ANOVA model.

Landsat and Phenocam Relationship
Due to the spatial resolution (30 × 30 m) of the Landsat data, metrics were calculated and used at the whole meadow level, instead of at individual community levels. Across all meadows, treatments, and years, Pearson's r values demonstrated a strong correlation between Landsat NDVI and phenocam GCC, and RMSE values demonstrated strong similarities in the prediction errors of the comparisons (Table 5). RMSE values never rose above 0.108, and Pearson's r values never dropped below 0.91. For both 2019 and 2020, the meadow with the highest Pearson's r value was Lower B (0.983 and 0.984). The meadow with the lowest Pearson's r value for both years was Upper A (0.911 and 0.912).

Soil Moisture
Determining the RD dates for the soil moisture of the meadows demonstrated the effect year-to-year differences had on meadow soil moisture retention (Figure 9). The average RD DOY for 2019 was 190, while the average DOY for 2020 was 144. The treatment with the latest DOY in 2019 was the grazed (192 average), and the treatment with the latest DOY in 2020 was the uncontrolled (162).

Influences of Climatic and Community Variables
Prism data showed that the study encompassed both the wettest year over the past decade in the meadow complex (2019, 419 mm annual precipitation) and the driest (2020, 175 mm annual precipitation). Interestingly, the effect this change in precipitation regime had on plant phenology was more distinct in the strictly camera-based data than in the random forest (RF) model, which combined both field-and camera-based observations. From the RF model, year to year variability was designated as the second least important variable (as determined by its mean decrease in Gini). In contrast, soil moisture differences

Influences of Climatic and Community Variables
Prism data showed that the study encompassed both the wettest year over the past decade in the meadow complex (2019, 419 mm annual precipitation) and the driest (2020, 175 mm annual precipitation). Interestingly, the effect this change in precipitation regime had on plant phenology was more distinct in the strictly camera-based data than in the random forest (RF) model, which combined both field-and camera-based observations. From the RF model, year to year variability was designated as the second least important variable (as determined by its mean decrease in Gini). In contrast, soil moisture differences within the season of growth had one of the highest levels of importance. When analyzing year to year differences picked up by camera-derived phenophases, the changes were more perceptible, with the dry year of 2020 having consistently earlier phenophases across the meadows, which is expected due to earlier snow melt. It should also be noted that within the RF model, temperature variability within the growing season was much less important in determining phenological stages than DOY, life form, and soil moisture.
The type of vegetation community also played a less important role in determining phenological stage than expected within both the RF model and the ANOVA model for the field data. The ANOVA model only showed community significance for the S phenology stage, and the RF model had community as the sixth most important variable (out of 11). While delineation of these community types was determined to be fairly accurate, overlap may have occurred leading to error. Encroachment of species from one community type into another was common, and is clearly demonstrated in Figure 4 [55,56]. Conversely, when comparing camera-derived dates, the ANOVA showed significant differences between communities for all phenophases except UD. The major differences in phenophases were seen between mesic vegetation and the other two community types. For both years, very little variation was seen in the UD dates (probably due to high soil moisture in all meadows early in the year). In 2020, however, as the growing season progressed, the phenophase dates of the mesic communities came earlier. This could have been due to the mesic community soils drying out quicker than the wet community soils, as well as the woody shrubs in the drier communities having deeper roots and naturally later seasons of growth than the mesic sedge communities. Regardless, for the earlier portion of the year, there was not much difference between vegetation communities. It should also be noted, that while phenophase UD dates started at roughly the same time for both years across all vegetation community types, the greatest disparity between years was at the downturn and end of growing season (i.e., DD and RD) in the mesic communities (Figure 7). This demonstrates that the mesic communities were not only affected the most by interannual environmental changes, but also by year-to-year environmental changes. The clarity provided by the phenocams shows the importance that multiple temporal resolutions of observation can have on understanding phenology in these complex ecosystems.
The RF model as a whole, had issues with differentiating between phenological stages within the meadow complex. While the model was fairly accurate in predicting when vegetation reached stages L and D, the error linked with stages I and S was significantly higher. This may largely be explained by the difficulty associated with properly identifying the two stages in the field. Certain species required greater observer experience in recognizing phenological shifts and human error may be present in the data. Dense sagebrush cover may have also obscured vision and made it difficult to ascertain field metrics of understory forbs and graminoids. Additionally, variability in genotypes and small-scale variability caused by competition or edaphic factors could lead to differences of phenology of individual species [57,58]. This could have allowed for a high level of irregularity in individual species' flowering and seeding time, irregularity which the model may not have detected.
It should be noted, as well, that when the two stages were differentiated, they usually occurred fairly close together in the context of the growing season as a whole. The dates where the highest percentage of observed species were in the I and S stages were very close across all meadows, while the dates for L and D were much further from the next closest stage. This would make it very difficult for the model to separate the I and S stages from one another. In future iterations of the model, the I and S stages could possibly be combined into a single group. While more basic, this would still help managers and researchers in understanding when plants would be susceptible to grazing and could influence management decisions.

Grazing Impact on Phenology
One of the major limiting factors of this study was the deficient number of meadows that had an uncontrolled grazing treatment. This was due, in large part, to time and monetary restrictions. When analyzing the effects grazing had on phenology, the RF model seemed to lack the ability to pick up variation across treatments. This was demonstrated by grazing having the third lowest level of importance among the variables used. Cameraderived phenophases and field-observed phenology dates also showed little variation between grazing treatments. Some variation was seen within the RD dates, but this was not deemed significant. With more years of grazing management, and more meadows included in the uncontrolled treatment, differences could potentially be seen.
We theorize that any differences seen between the grazing treatments in the future may be partially explained by soil moisture being retained in the soil later in the season. We fitted double logistic curves comparable to the ones fitted for GCC to calculate the "phenophase dates" for the soil moisture values in each of the meadows (Figure 9). The uncontrolled meadow had a similar RD DOY for soil moisture in 2019 to the other meadows, but in 2020 the RD DOY shifted much later into the year when compared to the others. Zhang et al. [59] hypothesized and observed that increased grazing in alpine meadows initially led to lower levels of transpiration from vegetation, allowing for moisture to stay later in the soil. Eventually, however, long term grazing led to an increase of bare ground and evaporation, limiting soil moisture in comparison to less grazed areas. Similar trends were shown in other studies across rangeland ecosystems [60][61][62]. If similar grazing patterns were to continue in Haypress, a shift in soil moisture to earlier RD dates could potentially occur within the uncontrolled meadow. Alternatively, sustained grazing has the potential to delay plant development (i.e., reduce plant size), and in so doing limit water uptake until later in the season. In a study conducted by Virgona et al. [63], they demonstrated that prolonged grazing in wheat fields could potentially elongate the growing season and delay development. They stated that within their experiment, for every 4-5 days grazed, antithesis and maturity was delayed by 1 day. This and other studies demonstrate how grazing can affect the maturation of vegetation [64][65][66]. This delayed necessity for water uptake could be another explanation for why soil moisture in the uncontrolled meadow was higher when compared to other meadows and may continue to be in future years.

Relationships between Spatial Resolutions
Phenocam derived GCC demonstrated a strong correlation with field observations by being the variable with the second highest degree of importance in the RF model. Additionally, Landsat NDVI maintained a high level of correlation with phenocam GCC across all years, treatments, vegetation communities, and sites observed (see supplementary Figure S1 for more detail). These strong relationships make a clear argument for the ability to develop statistical models which could extrapolate phenology metrics to other Great Basin Meadow systems using Landsat. With a wider sample size of meadows, years, and grazing treatments, and more effort to define and model the relationship between phenocam and satellite imagery, Landsat data could also potentially be used to extrapolate grazing pressure to other meadow systems.

Conclusions
Overall, this study shows how remote sensing and on-the-ground science can work in coordination to create a deeper understanding of phenological changes across gradients. Within the Haypress complex, more years of observation need to be incorporated in order to examine vegetational shifts in years that do not embody opposite ends of its climatic spectrum. Furthermore, additional meadows could be included to analyze how heavily grazed systems change in relation to properly managed environments and see if these systems emulate similar patterns to previous studies. Understanding the influence grazing has on these areas, in the context of phenology, will allow for the implementation of management strategies that do not create dissonance between the timing of meadow resources and the needs of those that rely on them.