Ecosystem Gross Primary Productivity After Autumn Snowfall and Melt Events in a Mountain Meadow

Vegetation productivity is increasing in the U.S. Northern Great Plains but decreasing in some nearby Northern Rocky Mountain grasslands due to increases in aridity. It is unclear if decreases to montane grassland productivity from drying autumns can be partly offset by late‐season green‐ups after precipitation events. These include the multiple snowfall and snowmelt periods that often characterize the summer‐to‐winter transition, but are difficult to observe due to logistical constraints. Here, we quantify changes to vegetation indices and ecosystem carbon uptake that occurs after snowfall and melt in climatological autumn in a montane grassland in Montana, USA using eddy covariance, phenological camera, and remote sensing analyses. Carbon dioxide flux follows a diurnal pattern after autumn snowmelt events despite overall ecosystem C loss, suggesting that post‐snowmelt photosynthesis helps dampen C loss during autumn and provides fresh photosynthate to support ecosystem functioning. Light‐saturated photosynthesis after three snow events was not statistically different than before snowfall (∼6 μmol CO2 m−2 s−1 in 2016 and ∼2.5 μmol CO2 m−2 s−1 in 2017). Observations are consistent with the notion that canopy photosynthesis is resistant, rather than resilient, to autumn snow. MODIS indicates that post‐snow greenups can occur, but not every year. Such events likely play a small role in the annual ecosystem carbon balance but may be disproportionately important for organisms faced with dwindling late‐season forage. Future efforts should seek to understand the community and ecosystem consequences of vegetation functioning during autumn as part of an expanded effort to understand phenological changes during this under‐studied and changing time of year.

Montane grasslands support diverse biological communities (Körner & Spehn, 2019) and provide critical ecosystem services (Körner, 2004) that are under threat from changes to land use and climate (Egan & Price, 2017) including snow dynamics (Harpold et al., 2012). Carbon sequestration is a critical ecosystem service that may be substantial in montane grasslands (Reed et al., 2021) but is highly variable across space and time-due in part to climate and management pressures (Rogiers et al., 2004;Schmitt et al., 2010;Tenhunen et al., 2009;Wohlfahrt, Anderson-Dunn, et al., 2008)-which may result in net annual CO 2 losses (Rogiers et al., 2008). Regardless of net annual carbon efflux, the documented pronounced late season declines of vegetation productivity in montane grasslands of the northern Rocky Mountains is cause for concern for the organisms that rely on fresh vegetation for overwintering success (Hurley et al., 2014(Hurley et al., , 2017Loe et al., 2020).
Despite large declines in biomass during the late growing season before snowfall in a mountain meadow in Montana , research at the site and communication with land managers has led to the speculation that precipitation pulses-often from snowfall and subsequent ablation (hereafter "melt")-support a "recovery" of production during the transition from the dry summer season to winter. This notion is consistent with a hypothesis that montane grassland photosynthesis is resilient to late-season snow events and can quickly re-establish previous function after melt, much in the same way that grass photosynthesis tends to re-establish quickly following harvesting (Novick et al., 2004;Stoy et al., 2008), grazing (Parsons et al., 1983), and heat waves (Hoover et al., 2014). Alternately, canopy photosynthesis in mountain meadows could simply be resistant to snow disturbances following a growing body of cold-season research that demonstrates that photosynthesis can occur under snow (Lundell et al., 2008;Starr & Oberbauer, 2003;Street et al., 2012) and the notion that grass productivity is resistant to minor disturbances like mild drought (Batbaatar et al., 2021). In this scenario, ecosystem carbon uptake would not meaningfully differ before and after snow events and ecosystem carbon uptake would have little response to snow apart from the limitations to photosynthesis caused by cold conditions and diminished light availability under the snowpack. Most research to date has characterized the rapid re-establishment of photosynthesis after snowmelt in spring (Julitta et al., 2014;Vitasse et al., 2017;Wohlfahrt, Anderson-Dunn, et al., 2008), with fewer studies of ecosystem carbon uptake in response to short-duration snow events during the autumn-to-winter transition.
Ecological research in montane ecosystems during autumn and the transition from the vegetative growing season to winter is lacking in part due to the difficult travel conditions created by snowfall and the subsequent melt periods that make roads impassable and over-snow transport onerous. This has hindered our understanding of ecosystem functioning during this critical and changing period (Piao et al., 2008), especially compared to its spring counterpart (Gallinat et al., 2015). To help remedy this situation for the case of montane grasslands, we used automated systems to quantify the effects of autumn snowfall and melt periods on ecosystem carbon uptake and loss in a montane grassland in Montana, USA that has been studied for decades. Specifically, we characterize the surface-atmosphere exchange of carbon dioxide and energy, greenness, and the normalized difference vegetation index (NDVI) before, during, and after snow accumulation and melt events during climatological autumn (September, October, and November) using eddy covariance, phenocam, and remote sensing measurements and ask if observations follow a pattern of ecological resilience or resistance to autumnal snow.

Site Information
Observations were made at the Bangtail Study Area Weaver & Collins, 1977), which hosts the Bangtail Mountain Meadow eddy covariance research site (US-BMM, Stoy, 2021). The study site is located 24 km northeast of Bozeman, MT, USA on a windswept ridge at 2324 m above sea level at 45.7830 N, 110.7776 W. The site served as the high elevation grassland biome site for the International Biological Program (1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974) and hosts the longest-active snow manipulation experiment of which we are aware Yano et al., 2015). The mean annual temperature is 2.3 °C and the mean annual precipitation is 981 mm as estimated by daymet (Thornton et al., 2014). Vegetation is dominated by Idaho fescue (Festuca idahoensis) with increasing cover of the invasive Kentucky bluegrass (Poa pratensis) and soils are loamy-typic cryoborolls (Weaver, 1977;Weaver & Collins, 1977;Yano et al., 2015). The study area has been fenced against cattle since the 1930s but is available for forage by wild large ungulates (e.g., mule deer (Odocoileus hemionus) and elk (Cervus canadensis)) and smaller subalpine animals including marmot (Marmota caligata), deer mouse (Peromyscus maniculatus), northern pocket gopher (Thomomys talpoides), and more.

Eddy Covariance and Micrometeorological Instrumentation
A 3-m tall tripod tower was installed on 30 September 2016 to host eddy covariance and micrometeorological instrumentation ( Figure S1 in Supporting Information S1). Measurements became available on 1 October 2016. Eddy fluxes were measured using a three-dimensional sonic anemometer (CSAT-3, Campbell Scientific, Logan, UT, USA) mounted at 153 cm above ground level and an enclosed infrared gas analyzer (EC-155 Campbell Scientific)-both recording at 10 Hz to a CR6 datalogger (Campbell Scientific)-as part of a CPEC300 eddy covariance system. Air temperature (TA) and relative humidity were measured using an HMP45C (Vaisala, Vantaa, Finland) at 200 cm and used to calculate the atmospheric vapor pressure deficit (VPD). Incident and outgoing short-and longwave radiation were measured using an NR01 net radiometer (Hukseflux, Delft, The Netherlands) at 176 cm. Outgoing longwave radiation observations were used to estimate the radiometric surface temperature (T surf ) using the Stefan-Boltzmann equation, assuming a view factor of unity and an emissivity of 0.98 to account for all surfaces observed by the radiometer including grass, bare ground, and snow (Snyder et al., 1998). The shortwave albedo (ALB) was calculated as the ratio of outgoing to incident shortwave radiation, and we use the median of measurements taken between 1100 and 1300 local time to minimize solar zenith angle effects. The normalized difference vegetation index (NDVI) and photochemical reflectance index (PRI) were measured at 176 cm using Spectral Reflectance Sensors (Decagon (now METER), Pullman, WA, USA). We likewise use the median of measurements taken between 1100 and 1300 to calculate daily values. Soil heat flux (G) was measured using two HFP01SC heat flux plates (Hukseflux) at 10 cm below the soil surface. Three CS650 soil temperature (T soil ) and moisture content (SWC) sensors (Campbell Scientific) were positioned horizontally at 10, 20, and 40 cm below the soil surface to account for a rooting depth observed to be at least 50 cm. An additional CS650 was positioned vertically to make an integrated measurement of the top 30 cm of soil. We present measurements from the sensor at 10 cm to characterize conditions near the surface including potential freezing. An SR50 sonic distance sensor (Campbell Scientific) was pointed downward at 146 cm to measure distance to surface; we use the median observation from each day to characterize vegetation height or snow depth (h). Micrometeorological observations were measured every second, averaged to half-hour periods, and saved to a CR1000 datalogger (Campbell Scientific). A list of abbreviations is presented in Table 1.

Eddy Covariance Data Processing
Half-hourly measurements of the eddy flux of carbon dioxide (Fc), latent heat flux (LE), sensible heat flux (H), and momentum flux were calculated using EddyPro (LiCor, Lincoln, Nebraska, USA). To calculate fluxes, we used double rotation (Kaimal & Finnigan, 1994) rather than the planar fit method due to rapid changes in canopy and snow height, as well as block averaging, covariance maximization for time lag detection, and highand low-pass filters (Moncrieff et al., 1997(Moncrieff et al., , 2005. Spike removal was performed as described by Vickers and Mahrt (1997) and defined as more than 3.5 standard deviations from the mean mixing ratio for water vapor and carbon dioxide and five standard deviations for vertical wind velocity. Fluxes that did not pass the Mauder and Foken (2011) quality control tests were excluded from the analysis. We assumed that storage of carbon dioxide below the 153 cm instrument height was negligible such that Fc approximates the net ecosystem exchange of carbon dioxide (NEE). Open access data from US-BMM are published on the Ameriflux website at https:// ameriflux.lbl.gov/sites/siteinfo/US-BMM.
NEE was partitioned into its components, gross primary productivity (GPP) and ecosystem respiration (RE), using REddyProc (Wutzler et al., 2018). We explored flux partitioning methods that use nighttime NEE measurements (Reichstein et al., 2005) and both nighttime and daytime NEE measurements (Lasslop et al., 2010) given that results from these two methods can differ and impact the interpretation of GPP (Reichstein et al., 2012;Stoy et al., 2006). The partitioning method based on nighttime data develops a relationship between RE and air temperature as described in Reichstein et al. (2005): where R Ref (μmol CO 2 m −2 s −1 ) is the reference respiration at a reference temperature, T Ref , set to 15°C and T 0 is set to −46.02°C (Lloyd & Taylor, 1994;Reichstein et al., 2005;Wutzler et al., 2018). Nighttime data are constrained between local sunrise and sunset and when incoming shortwave radiation is less than 10 W m −2 . E 0 is estimated using 15-day windows and aggregated to an annual estimate, and a temporally varying R Ref is estimated using the annually aggregated E 0 (Wutzler et al., 2018). The relationship developed from nighttime data is extrapolated using TA measurements during the day to estimate daytime RE. GPP is then calculated as: GPP = RE − NEE.
The partitioning method based on nighttime and daytime data estimates NEE using a light response curve following Lasslop et al. (2010): where α NEE (μmol CO 2 J −1 ) is the canopy light use efficiency before light saturation is reached, β NEE (μmol CO 2 m −2 s −1 ) is the maximum CO 2 uptake rate at light saturation, R g is the incoming shortwave radiation at the surface (W m −2 ), and γ (μmol CO 2 m −2 s −1 ) is RE at a R g of 0 W m −2 β NEE in the Lasslop et al. (2010) algorithm is calculated as: where β 0 is the maximum CO 2 uptake rate at light saturation without VPD constraints and β NEE = β 0 when VPD ≤ VPD 0 . VPD 0 is set as 10 hPa. E 0 is estimated using nighttime data in moving windows, and the resulting estimates are smoothed (Wutzler et al., 2018). Prior estimates for R Ref are also estimated using nighttime data for the moving window with smoothed E 0 . Finally, the parameters (R Ref , α NEE , β 0 , k) are estimated using daytime data for each time window (Wutzler et al., 2018).
The first term on the right-hand side of Equation 2 is related to the rectangular hyperbola that can be used to represent the response of GPP to PPFD and we estimated its parameters using nonlinear least squares (nls) in R (R Core Team, 2021) to obtain β for periods before and between snow events, called "inter-snow periods" (IS). We used ANOVA with the post-hoc Tukey's HSD test to calculate differences in β among different measurement periods at the 0.05 significance level, using degrees of freedom to account for the number of observations and applying the Bonferroni correction for multiple comparisons. The calculation of the 95% confidence intervals was likewise adjusted for multiple comparisons following Dunn (1961).
The β parameter of the rectangular hyperbola (Equation 4) can reach large, unobserved values of GPP at large, unobserved values of PPFD which may poorly describe ecosystem CO 2 flux responses to light (Reichstein et al., 2012). To ensure that our results are not sensitive to choice of light response model, we repeated the analysis using the Mitscherlich model ) where the initial slope of the light response curve and GPP at saturation have similar meanings to their counterparts in Equation 4 but will take different values and are distinguished with subscript "M." Latent heat flux (LE) and sensible heat flux (H) were also gapfilled using REddyProc, which applies mean diurnal sampling to estimate missing values of these quantities (Falge et al., 2001).

Phenocam Observations and Analysis
A 5 MP NetCam SC phenocam (StarDot Technologies, Buena Park, CA) was installed 75 m SW of the eddy covariance tower during August 2015. The phenocam pointed north and captured an area that is typically upwind of the tower given the predominantly westerly winds (Figure 1). Observations were available for November and December 2016 after eddy covariance observations began, and for 2017 until early October. and For this analysis, we use the daily 90th percentile values of the GCC and RCC as recommended by Richardson et al. (2018) to minimize the effects of diurnal variability in illumination geometry and its effects on phenocam imagery.
Phenocam observations were infrequently available when flux measurements were made, and we focus our phenocam analyses on the 2017 autumn period.

Overview
We first present micrometeorological, surface-atmosphere exchange, and radiometric data for the 2016 study period, which we define as October 1 when measurements began until November 12 when flux measurements became unavailable due to power system failure and challenging road conditions made it difficult to repair. We follow the 2016 results with the 2017 findings and define a study period of September 10-four days before first snowfall for consistency with the 2016 study period-until November 12. Measurements were available until 12 December 2017, but with consistent snow cover after October 29 as we detail later. We then discuss measurements made during the summer of 2017 for comparison with autumn measurements and describe MODIS observations for the years prior to, during, and after the eddy covariance study period.

Observations During Autumn, 2016
G decreased from a maximum of 34 W m −2 when measurements began on 1 October 2016 ( Figure 2) and did not exhibit a diurnal cycle 4 days later, when T surf was equal to or less than 0 °C, SWC increased from 30% to 35%, ALB increased to values greater than 50%, NDVI decreased to less than 0, and h (here suggesting changes in snow height) increased from 0.12 to 0.20 m (Figure 3). These observations are consistent with snow presence, and we subsequently use micrometeorological, radiometric, and flux observations to determine snow presence or absence (Table 2).
NEE exhibited a distinct diurnal signal before the first snowfall event of autumn, 2016, that was dampened when snow was present (Figure 4). The nighttime partitioning approach (Reichstein et al., 2005) estimated that GPP was usually no greater than 2 μmol m −2 s −1 during peak daytime periods when snow was present. The daytime partitioning approach (Lasslop et al., 2010) frequently estimated a GPP of ∼2 μmol m −2 s −1 during peak daytime periods when snow was present ( Figure S2 in Supporting Infomation S1). From these observations, nighttime partitioning is arguably more conservative for estimating GPP in the presence of snow. We focus the rest of our analysis on GPP and RE values obtained using the Reichstein et al. (2005) nighttime partitioning method.
We identified four snow periods in 2016 (Table 2, Figures 2-4) before measurements became unavailable for the rest of the year on November 12. After the first snowfall, T soil was lower and SWC higher than the pre-snow period for the remainder of the 2016 study period. T soil remained lower than 6 °C after first snowfall for the remainder of the study period, but TA consistently reached 10 °C or greater and T surf was frequently greater than 20 °C during midday periods between snowfall events. Mid-day reflectance measurements (NDVI and PRI) did not differ for the periods before and after snow (Figure 3). Mid-day GPP decreased from ∼6 μmol m −2 s −1 before the first snowfall to between ∼2 and 4 μmol m −2 s −1 after the fourth snow event (Figure 4). The statistical analysis of light response curve parameters indicated no significant difference in β for the period before snow and first three inter-snow periods, which averaged 6.14 μmol CO 2 m −2 s −1 (Figure 5, Tables 2 and S1 in Supporting Information S1). β during the third inter-snow period was significantly lower (4.62 μmol CO 2 m −2 s −1 ) than the first and second intersnow periods and the fourth inter-snow period was significantly lower than all other snow-free periods (3.44 μmol CO 2 m −2 s −1 ) ( Figure 5, Tables 2 and S1 in Supporting Information S1). Similar results were obtained using the Mitscherlich model with the exception that β calculated for the third intersnow period was lower than before snow and the first two inter-snow periods (Table S3 in Supporting Information S1). RE decreased from 3 μmol m −2 s −1 to 1 μmol m −2 s −1 throughout the 2016 autumn measurement period ( Figure 4). As a consequence, the ecosystem was a net sink of CO 2 (i.e., NEE was negative) during midday of all snow-free periods during the measurement record in 2016 (Figure 4) but was a source of C to the atmosphere during other times of day and lost ∼50 g C m −2 from October 2 until 12 November 2016.

Observations During Autumn, 2017
Autumn 2017 was characterized by two weeks of snow presence from September 15 until September 29, followed a three-day snow-free period, then four days of snow, then nearly a week without snow, and additional snow and no-snow periods of less than week in duration before consistent snow presence began on October 29 (Table 3). We determined these dates using a combination of G, T surf , turbulent flux observations, and phenocam observations when available, and we discuss challenges in determining snow presence in the Discussion.
SWC was ∼10% before snowfall in 2017 ( Figure 6)-recalling that SWC was ∼30% before snow fell in 2016 (Figure 2)-and increased to 35% after  Note. Light response curves were not calculated for snow periods despite non-zero GPP (e.g., Figure 4) due to the frequent patchy snow conditions (e.g., Figure 1). *,**,*** indicate statistical significance at the 0.05 significance level after applying the Bonferroni correction Additional statistics are provided in Table S1 in Supporting Information S1. a Measurements became unavailable for the rest of the year after 12 November 2016.

Table 2 The Dates of the Before-Snow Period, Snow Presence (S) and Inter-Snow Periods (IS), Numbered Sequentially, at the Bangtail Mountain Meadow (US-BMM) Eddy Covariance Research Site in Montana, USA During 2016 With Degrees of Freedom (df), Light-Saturated GPP (β), and Its Standard Error
the first 2017 snow event. TA intermittently reached nearly 20 °C during the periods between snow in 2017, as did T surf , and T soil remained above 0°C until November 12 when the study period ended for consistency with data availability during 2016 ( Figure 6). As opposed to 2016, NDVI stayed consistently at or below 0 and ALB remained above 50% during the periods between snow presence (Figure 7). The phenocam was active during the first and second snowfall and melt periods in 2017, and RCC declined sharply from nearly 0.6 to less than 0.35 after the first snow event (Figure 8). RCC recovered to values approaching 0.4 during the first inter-snow period and nearly 0.45 after the second snow period after which measurements became unavailable. GCC was less sensitive to snow presence and absence (Figure 8).
Maximum daily GPP before and after the first two snow events was on the order of 3 μmol m −2 s −1 (Figure 9); recall that light-saturated GPP (i.e., β) reached values of ∼6 μmol m −2 s −1 in 2016 (Figures 4 and 5). β was 2.47 μmol m −2 s −1 before snowfall and was not significantly different before snow or during the intersnow periods but decreased to 1.50 μmol m −2 s −1 during the third inter-snow period ( Figure 10, Tables 3 and S2 in Supporting Infomation S1). Similar results were found when using the Mitscherlich model (Table S4 in Supporting Information S1). As opposed to 2016, NEE was rarely less than 0 μmol m −2 s −1 such that the ecosystem was almost continuously losing carbon after the first snowfall in 2017 (Figures 4 and 9). The ecosystem was a consistent net source of C to the atmosphere during climatological autumn as a consequence. The ecosystem lost 51 g C m −2 during the period from October 1 to 12 November 2017, similar to 2016. There was a notable increase in NEE and RE from values less than 2 μmol m −2 s −1 to values greater than 2 μmol m −2 s −1 around 17 September 2017, a snow-covered period when TA was greater than 0°C, SWC increased from 10% to 20%, LE reached values of nearly 100 W m −2 , and H was negative during the day indicating heat flux into the snowpack (Figures 6 and 9).

Observations During Summer, 2017
Eddy covariance and micrometeorological measurements were available beginning in mid-July 2017 onward (Figure 11), and we present these observations to place the 2017 autumn measurements in context. Light-saturated GPP (i.e., β) was 18.6 μmol m −2 s −1 during late July (Figure 12), already decreased to 9.02 μmol m −2 s −1 as drought progressed in early August, decreased further to 6.16 μmol m −2 s −1 during late August, and reached  ∼3 μmol m −2 s −1 during early September as noted above ( Figure 10, Table 3). In other words, β was approximately 1/6th of its maximum observed value in summer before snowfall in 2017.

MODIS
We analyzed NDVI MODIS observations greater than 0.4 to help avoid noisy measurements at lower values that are consistent with full or patchy snow cover. A notable increase (i.e., greening) in NDVI MODIS occurred in 2015 from 0.45 on October 7 to 0.5 on October 27, and NDVI MODIS approached 0.5 in 2017 after a long period of low NDVI MODIS values with an inconsistent vegetation signal that began on September 12 ( Figure 13). Other years, namely 2016 and 2019, did not demonstrate any late season increase in NDVI MODIS although measurements in 2016 rarely dropped below 0.4 during the latter part of the calendar year.

Overview
Decades of research at the study site and communication with local ranchers led to the suggestion that canopy photosynthesis would "recover" after drought-breaking autumn precipitation, which often falls as snow. The ecosystem was not likely under drought stress before snow fell in 2016 as SWC was ∼30% (Figure 3) but was only ∼10% before snowfall during 2017 (Figure 7), well below the point at which grasses are often considered drought stressed (Rodriguez-Iturbe et al., 1999;Rodríguez-Iturbe & Porporato, 2007), and consistent with declining GPP (Figures 9, 11, and 12). Observations from 2017 indicate that GPP was rarely greater during inter-snow periods than it was before snowfall despite added moisture (Figures 4 and 9). β was not significantly different before or after the first three snow events in either year (Figures 5 and 10, Table 2), following the statistical analyses that included Bonferroni corrections to significance values and confidence intervals (Dunn, 1961). Observations do not necessarily suggest that there is a lag in the recovery of GPP after snow melt nor notable regrowth (Figures 4 and 9). In other words, rather than a recovery of photosynthetic function and thereby ecological resilience, our observations suggest that the vegetation probably never lost photosynthetic

Note.
A light response curve was not calculated for the single day of available data during IS5, and the data were not significantly related to the light response curve during IS4, which is excluded from the statistical analyses. * indicates statistical significance at the 0.05 significance level after applying the Bonferroni correction Additional statistics are provided in Table  S2 in Supporting Information S1. a Measurements became available on 12 July 2017, but for consistency with the 2016 study period we fit light response curves using data collected for the four days preceding the first snow event beginning on 10 September 2017. b Measurements became unavailable for the rest of the year after 18 December 2017.

Table 3
The Same as Table 2 But for 2017 Figure 6. The same as Figure 2 but for 2017. Periods identified as having snow (Table 3)  function under snow in the first place, at least after the first three snow events. Instead, observations are consistent with the notion that autumn snowfall is a minor disturbance against which mountain meadow photosynthesis is resistant.
Temperate, boreal, and arctic grasses can photosynthesize at temperatures at or below 0°C (Skinner, 2007;Tieszen, 1973;Tieszen & Sigurdson, 1973), and some grasses maintain photosynthesis under cold and dark conditions that simulate the subnivean environment (Höglind et al., 2011). Our observations suggest finite GPP even in the presence of snow (e.g., Figures 4 and 9). At the same time, the study site has diverse vegetation, and it is unclear if grasses are the primary contributor to autumn GPP. Some mosses can photosynthesize under snow (Street et al., 2012) and are characterized by light-saturated photosynthesis (i.e., β) ranging from nearly 0 to 7 μmol m −2 s −1 (Skre & Oechel, 1981;Street et al., 2012;Van Gaalen, Flanagan, & Peddle, 2007), similar to our observations during periods between snow events in 2016 ( Figure 5). At the same time, moss cover at the site is sparse and likely insufficient to explain observed β values on the order of 7 μmol m −2 s −1 on a ground-area basis. Instead, photosynthesizing vascular vegetation-that often occurs toward the bottom of grass canopies in autumn (Migliavacca et al., 2011)-is more consistent with available observations (see Figure S1 in Supporting Infomation S1). These patterns are often difficult to discern with phenocams that point laterally (Figures 1 and 8) (Migliavacca et al., 2011). Tower-mounted NDVI sensors point downward, and measurements often exceeded 0.2 during periods between snow (e.g., Figure 3). These are low values for grass NDVI (Gamon et al., 1995) that are consistent with the small carbon fluxes that we observe (Figures 4,5,9,and 10). Carbon uptake is also lower at our montane study site than most global grasslands; peak light-saturated NEE at 1500 μmol m −2 s −1 PPFD following the approach of (Ruimy et al., 1995) is 10.9 μmol m −2 s −1 at our study site, versus a global average of 25.5 μmol m −2 s −1 (Gilmanov et al., 2010).
Measurements gave no evidence that mountain meadow GPP was greater after snow than before snow, even when more moisture and favorable growing conditions occurred after melt. It was unclear if new leaves flushed, but no observable height growth occurred between snow events (Figures 3  and 7). We cannot say for certain which grass species were responsible for photosynthesis after snow, but observations in nearby valleys suggest that the invasive Poa pratensis rather than native Festuca idahoensis is more likely to remain green in the cold season, with consequences for ecosystem function   (Table 3).
as invasive grasses continue to establish at the study site. A comprehensive monitoring framework would require species-level studies to help reconcile population and community dynamics that we were unable to observe with certainty. Species-level measurements after snow at the study site require travel during uncertain conditions, and/ or high-resolution imagery that is increasingly used in conjunction with machine learning approaches to estimate grass growth and phenology (Oliveira et al., 2020;Viljanen et al., 2018).

Ecosystem Dynamics During Snow Presence
A notable event occurred around 17 September 2017, when snow was present, but TA reached nearly 10 °C and H became negative, indicating heat flux into the snowpack (Stoy et al., 2018), but T surf remained at or near 0 °C indicating that snow did not fully melt. SWC increased from 10% to 20% during this event and NEE also increased, suggesting enhanced soil respiratory activity. We posit that the snowpack was melting during this period but did not fully ablate, and the input of soil water enhanced RE owing to the Birch effect (Jarvis et al., 2007): the rapid increase in decomposition rate and biogeochemical activity induced by rewetting dry soils. We further speculate that the soil water input and RE were synchronized in time because of rapid pressure pumping (i.e., enhanced diffusion) through snow at this characteristically windy montane site (Berryman et al., 2018;Bowling & Massman, 2011;Rains et al., 2016). Decadal snow fence observations at the site suggest that enhanced snowpack duration depletes soil C, N, and P (Weaver & Collins, 1977;Yano et al., 2015), but the microbial mechanisms of potential under-snow Birch effects in ecosystems that receive drought-breaking snow need to be substantiated with direct observations of soil biogeochemistry.

The Challenges of Patchy Snow
We used micrometeorological and radiometric observations to identify snow presence using instruments on or near the tripod tower ( Figure S1 in Supporting Infomation S1), confirmed with phenocam measurements when available (Figures 1-3 and 6-8). Turbulent fluxes arise from the eddy covariance flux footprint, which was often in the phenocam scene that measured an area immediately west of the tower, the dominant wind direction (Figure 1). The physical tower infrastructure will cause snow deposition, and we suspect that this is the reason that 2017 melt events identified using micrometeorological and flux observations (Figures 6 and 9) differed from the radiometers (Figure 7). In other words, snow presence was often patchy (e.g., Figure 1)  (Table 3) are indicated in blue boxes. Figure 10. The same as Figure 5 but for 2017 with inter-snow periods (IS) as described in Table 3. and can be influenced by the measurement infrastructure itself by slowing wind speed to deposit snow. This leads to an observational challenge of identifying snow cover for the purpose of ecosystem-scale studies. Going forward, we propose a tripartite classification system with no snow, full snow, and patchy snow, which remains a challenge for observation systems (Kosmala et al., 2018). During melt, dynamic changes to the snow-covered area across space and time are being sampled by a dynamic flux footprint (Chu et al., 2021), and both can be highly uncertain. From our analysis such an observation system could not rely on tower-mounted instrumentation alone because the tower infrastructure can impact snow dynamics. In this case placing the tower within the field of view of the phenocam would be preferred.

Remote Sensing
Remote sensing is critical for understanding areas that are difficult to access like Northern Rockies grasslands (St. Peter et al., 2018). Differences are to be expected amongst RS sensors (Rossi et al., 2019) and different grassland sites , and the coarse spatial scale and diurnal return interval of MODIS gave different NDVI results than the tower-mounted instrumentation at our site (Figures 3, 7 and 13). Recent advances to geostationary satellite technology-for example, enhanced observational capacities in the visible and near infrared from the Advanced Baseline Imager on the Geostationary Operational Environmental Satellite (Khan et al., 2022)-provide new opportunities for understanding snow duration (Romanov et al., 2000;Romanov & Tarpley, 2004) and vegetation productivity after melt. Enhanced observational systems at multiple scales in space including tower-mounted instruments, phenocams, unpiloted aerial vehicles, and satellite remote sensing will further improve our understanding of autumn phenological changes in global ecosystems (Gallinat et al., 2015).

Conclusions
Grass photosynthesis tends to reestablish rapidly following harvesting (Novick et al., 2004;Stoy et al., 2008), grazing (Parsons et al., 1983), and heat waves (Hoover et al., 2014) if moisture is available, and the response is dependent on species and growth form . Grass photosynthesis also tends to establish rapidly after snowmelt in spring (Julitta et al., 2014;Vitasse et al., 2017), which we were not able to observe at our site due to equipment and power failure over winter. Our observations suggest that mountain meadow photosynthesis is resistant to autumn snowfall which appears to be a minor disturbance for snow events up to 2 weeks in duration. Understanding the resistance and resilience of grass vegetation to extreme climate conditions is  The relationship between GPP and PPFD for the periods July 12-31 ("late July"), August 1-14 ("early August"), August 15-31 ("late August"), and September 1-15 ("early September"), 2017, with fitted light response curves following Figure 5. critical as extreme drought (Hoover et al., 2014) and desertification (D'Odorico et al., 2013) threaten grasslands worldwide.
Snowpack dynamics are a key determinant of the phenology and productivity of subalpine meadows (Canaday & Fonda, 1974). We argue that spring snowmelt (or summer, depending on site) is only one of two periods of importance for montane grassland phenology and function. The transition to full snowpack, which may contain multiple snowfall and melt periods that can bring needed moisture to mountain meadows, is also important. Enhanced observational systems that consider population, community, and ecosystem-scale responses to ongoing climate changes during autumn will improve our understanding of the ecosystem services that montane grasslands provide across all seasons (Gallinat et al., 2015).

Acknowledgments
We acknowledge the Indigenous Nations on whose ancestral lands the study took place and recognize that the infrastructure used for this project was built on Indigenous land. We recognize multiple Indigenous Nations as past, present, and future caretakers of this land whose stewardship of the region was interrupted by their physical removal under the 1830 Indian Removal Act and through US Assimilation policies explicitly designed to eradicate Indigenous language and ways of being until the 1970s. We thank Dr. James Irvine for extensive technical support. Adam Cook, Dr. Geneva Chong, Justin Gay, Dr. Bill Kleindl, F. Aaron Rains, Dr. Angela Tang, and Bill Vandenberg provided additional technical support for which we are grateful, and we thank Dr. Susanne Wiesner for insightful comments on the manuscript.