AIC Akaike information criterion C carbon CWD climatic water deficit (mm) DBF deciduous broadleaf forest DoY day of year (1-366) ENF evergreen needleleaf forest GPP gross primary productivity (μmol m -2 s -1 ) GPP 1800 whole-forest photosynthetic capacity (GPP at high light) based on light response of NEE (μmol m -2 s -1 ) MAT mean annual temperature ( • C) MAP mean annual precipitation (mm) NEE net ecosystem exchange of CO 2 (μmol m -2 s -1 ) PAR photosynthetically active radiation (μmol m -2 s -1 ) R 2 coefficient of determination RMSE root mean square error (units vary) SWE snow water equivalent (mm) T air air temperature ( • C) VWC volumetric water content (m 3 m -3 ) seasonal metrics and time periods: SOS start of GPP season, the DoY in spring when GPP 1800 first reaches 10% of summer capacity (DoY) EOS end of GPP season, the DoY in autumn when GPP 1800 last reaches 10% of summer capacity (DoY) AS active season for GPP, the period between SOS and EOS, inclusive AS length length of active season for GPP -the number of days between SOS and EOS, inclusive spring ramp the period between SOS and 90% date in spring (the DoY when GPP 1800 first reaches 90% of summer capacity) autumn ramp the period between the 90% date in autumn (the DoY when GPP 1800 first reaches 90% of summer capacity) and EOS summer the period between 90% dates in spring and autumn, when forest is at peak photosynthetic capacity cumulative GPP during each season: ΣGPP AS sum of GPP during active season (g C m -2 ) ΣGPP spring sum of GPP during spring ramp (g C m -2 ) ΣGPP summer sum of GPP during summer (g C m -2 ) ΣGPP autumn sum of GPP during autumn ramp (g C m -2 ) SNODAS data products: SNODAS U.S. National Weather Service Snow Data Assimilation Program SAG the DoY when snow has fully disappeared (Snow All Gone, DoY) peak SWE annual maximum SWE (mm) day of peak SWE timing of annual maximum SWE (DoY) length of snowmelt difference between DoY of annual maximum SWE and SAG (number of days) PRLQ amount of liquid precipitation (rain, mm) PRSL amount of solid precipitation (snow, mm) solid fraction fraction of total annual precipitation falling as snow (%)
A complete understanding of how the terrestrial carbon (C) cycle responds to variability in environmental conditions is necessary for making accurate projections of the global C budget under future climate scenarios (Friedlingstein et al., 2022). Forecasting long-term C exchange of terrestrial ecosystems depends on understanding the environmental, biological, and biophysical controls of gross primary productivity (GPP), ecosystem respiration, and their balance (net ecosystem exchange of CO 2 , NEE) across seasonal and interannual timescales. This is a serious challenge, as the environmental controls of these fluxes are complex, and vary on a region-by-region or even case-by-case basis more often than universal relationships are found across ecological space (Baldocchi et al., 2018). In this study our focus is on GPP.
Across boreal and temperate evergreen and deciduous forests, controls of seasonal and interannual C exchange include phenological variability associated with the start of photosynthesis in spring (Richardson et al., 2009), fall senescence (Jeong et al., 2011;C. Wu et al., 2012a), or both (Desai et al., 2022;Goulden et al., 1996;Keenan and Richardson, 2015), environmental variability including seasonal temperature (Arain et al., 2022(Arain et al., , 2002;;Suni et al., 2003;Tanja et al., 2003), moisture (Goldstein et al., 2000;Thomas et al., 2009), light (Froelich et al., 2015), and disturbance (Aubinet et al., 2018;Finzi et al., 2020). In general, while carbon fluxes respond in sync with the environment on hourly to monthly time scales, the sensitivity of carbon fluxes to variation in weather progressively declines or becomes more difficult to detect seasonally and interannually (Richardson et al., 2007;Stoy et al., 2009;J. Wu et al., 2012), except in the case of extreme weather events (Zscheischler et al., 2014). It has been shown that less than half of interannual variability in NEE can be attributed to climatic factors, while the majority is due to variation in biological processes (Richardson et al., 2007;Shao et al., 2015). These considerations, combined with the effort and expense required to obtain multidecadal records, make understanding climatic influence on interannual C exchange challenging.
An important feature of ongoing climate change is reduced snow accumulation and related effects on water availability for plants. In western North America, snowpack reduction and earlier snowmelt are well documented (Hale et al., 2023;Mote et al., 2018;Siirila-Woodburn et al., 2021), and their decline is projected to continue (Barnett et al., 2005;Dierauer et al., 2019). Further, the amount of snow that melts intermittently during winter is increasing (Musselman et al., 2021). Continued reductions in snow may lead to longer seasons for photosynthesis. Here we refer to the photosynthetic period as the "active season" (defined formally in Section 2.3) and avoid the term "growing season", which is vague at best when considering the complexities of plant C allocation (Körner et al., 2023). Reduction in snow may also lead to drier soils and increased fire risk (Westerling, 2016). A recent study predicted that the number of snow-free days will increase from ~175 to ~250 by the end of the century in the central and northern Rocky Mountain region (Wieder et al., 2022), and warming and earlier spring snowmelt may extend the length of the active season. In many regions, snowmelt is a first order determinant of water availability (Barnett et al., 2005), and is particularly important for soil water infiltration and recharge of deep soil and groundwater (Jasechko et al., 2014). The timing of snowmelt tends to match the timing of peak soil water availability (Harpold and Molotch, 2015), which is important for transpiration (Cooper et al., 2020). Water from the winter snowpack may be used by plants well into the active season (Bailey et al., 2023;Goldsmith et al., 2022;Hu et al., 2010). Thus, the coupling of snowmelt and soil moisture is potentially important for interannual variability in C uptake and ecohydrological response to climate change.
Two alternate and competing ecological impacts of reduced snow amount and earlier melt have been proposed: the growth period effect and the moisture effect (Wang et al., 2018). A number of studies have investigated the implications of the growth period effect, which is defined as increased active season length due to earlier melt, and hence a longer period for growth. Some have reported increased C uptake with longer active seasons in temperate deciduous broadleaf forests (DBF, Goulden et al., 1996;Keenan et al., 2014;Richardson et al., 2009), as well as boreal DBF and evergreen needleleaf forests (ENF, Barr et al., 2002;Chen et al., 1999;Churkina et al., 2005). In DBF the potential for C uptake is constrained by new leaf emergence, therefore earlier start of the active season can increase the time the forest can function at maximum leaf-area. If soils are cold in the period after snowmelt, deciduous leaf emergence can be delayed (Desai et al., 2022), but warm soils during this period can offset potential C uptake benefits (Sanders-DeMott et al., 2020). In ENF, there is minimal seasonal change in leaf area, and instead photosynthetic function is subject to various environmental and biochemical constraints such as temperature and moisture availability, photoprotection, and photosynthetic downregulation, particularly in winter (Bowling et al., 2018;Chang et al., 2021;Monson et al., 2005;Verhoeven, 2014;Wolf et al., 2016). Therefore, in ENF, the ability to capitalize on earlier spring onset depends on whether additional constraints to photosynthesis are relieved.
While a longer active season can enhance GPP, increased spring or summer water stress due to reduced snow accumulation and/or earlier melt (the moisture effect) can outweigh the potential gains in GPP. A case study at a high-elevation subalpine forest near Niwot Ridge, Colorado (US-NR1) showed that earlier snowmelt leads to longer active seasons, but those years had decreased GPP due to late season soil moisture limitation (Hu et al., 2010). A more recent study at Niwot Ridge and other sites indicated that, in contrast to Hu et al. (2010), the interannual variability in net carbon exchange was not strongly related to active season length (Barnard et al., 2018). Stable isotope analyses have shown that in some regions, the water used by trees even late into the growing season originates primarily from soil moisture derived from snow melt (Allen et al., 2019;Berkelhammer et al., 2020;Hu et al., 2010;Martin et al., 2018;Phillips and Ehleringer, 1995). There is a growing consensus that increases in temperature and decreases in moisture associated with longer active seasons (the moisture effect) may decrease forest carbon sequestration (Knowles et al., 2018;Trujillo et al., 2012;Winchell et al., 2016).
While the moisture effect (also referred to as the seasonal compensation effect, Buermann et al., 2018) appears important in water-limited, snow-dominated ecosystems, such as the Rocky Mountains, water limitation is becoming increasingly characteristic of historically energy-limited temperate and boreal forests (Buermann et al., 2014;Butterfield et al., 2020;Denissen et al., 2022;Girardin et al., 2016;Peng et al., 2011). Terrestrial carbon cycle models overpredict the beneficial growth period effect and underpredict the adverse moisture effect that follows warmer springs (Buermann et al., 2018). Over the Northern Hemisphere, the strength and direction of the relationship between remotely-sensed snow and vegetation greenness is highly variable, and dependent on the relative dominance of the moisture and growth period effects (Wang et al., 2018). Additional remote sensing studies similarly show that the response of vegetation to changing snowpack is variable in magnitude and direction, and also spatially (Buermann et al., 2018;Xiong et al., 2019). Which of these controls dominates at any given site may be a function of ecosystem type, average moisture conditions, the legacy effect of snowmelt on summer soil moisture, seasonality of precipitation, and/or the degree of water limitation of vegetation.
Investigations of the impact of changes in snowpack on interannual carbon dynamics are unfortunately hampered by a lack of observations of snowpack characteristics at most flux towers. The water contained in the snowpack, referred to as snow water equivalent (SWE), is ecologically quite important. Combinations of ground-based observations that include SWE with remote sensing have led to progress in understanding interannual variation of forest greening (Knowles et al., 2017;Trujillo et al., 2012), but long-term SWE records co-located at eddy covariance flux tower sites are rare. Passive microwave remote sensing is quite helpful for estimating SWE (Kelly et al., 2003;Pulliainen et al., 2017), but at present has coarse spatial resolution, and accuracy is limited in the presence of forest canopies, deep snow, and mountainous terrain (Dozier et al., 2016;Mortimer et al., 2020;Vander Jagt et al., 2013).
As an alternative and/or addition to remote sensing, gridded climate reanalysis approaches can include assimilation of observational snow data and combine them with physical models to provide high-quality SWE estimates (Cho et al., 2020;Girotto et al., 2020;Zeng et al., 2018). This includes the Snow Data Assimilation System (SNODAS), developed by the US National Operational Hydrologic Remote Sensing Center (NOHRSC) and archived at the National Snow and Ice Data Center (NSIDC). SNODAS provides a 1 km 2 daily gridded estimate of SWE and related snow metrics over the contiguous US (since 2003) and southeast Canada (since 2010). SNODAS works by first ingesting data from the Rapid Update Cycle numerical weather prediction model, which are then downscaled and used to drive a physically based energy-balance and mass-balance snow accumulation and ablation model. The modeled output is then adjusted by data assimilation of all available ground, airborne, and satellite observations to produce a gridded estimate of daily SWE (Rutter et al., 2008). SNODAS generally works well to estimate SWE in forested areas (Artan et al., 2013;Clow et al., 2012), but has been primarily used for hydrologic applications. Its usefulness for ecological applications remains unexplored despite its relatively high spatiotemporal resolution.
In this study, we synthesized flux-tower observations of carbon fluxes and weather data from fourteen forest sites in the US and southeast Canada with gridded SWE and precipitation estimates from SNO-DAS. We used these data to study the potential legacy effects of snowpack dynamics on subsequent active season GPP, and which environmental controls might determine the relative dominance of the growth period vs. moisture effects of earlier snowmelt. We tested three hypotheses: H1) Earlier snowmelt leads to a longer active season for GPP, H2) Active season GPP is higher in years with longer active season, and H3) Active season GPP is dependent on peak SWE and timing of snowmelt (a winter to summer moisture legacy from the snowpack). These hypotheses provided a coherent framework to examine complex biophysical processes related to forest-atmosphere carbon exchange statistically, using linear regression and correlation analyses, and mixed effects models.
Very few flux towers include instruments to measure SWE or other snowpack parameters, so our analysis was limited to the region of SNODAS data availability (forests within the contiguous US and southeast Canada). Site selection criteria included seasonal snow cover, distinct periods of photosynthetic activity and dormancy, no recent disturbance, and 4 or more years overlap with SNODAS. Fourteen flux towers with 145 site-years of data met these criteria (Table 1). The forests are primarily in the Köppen-Geiger climate classification of Dfb (warm summer continental), with one exception classified as Dfc (US-NR1, subarctic/boreal, Peel et al., 2007). Mean annual air temperature ranged from 1.5 to 8 • C, mean annual precipitation 800-1250 mm, with mean annual maximum SWE (from SNODAS, Section 2.4) varying from 50 to 400 mm. The percentage of annual precipitation falling as snow (solid fraction, SNODAS) varied across sites 11-51% (Table B1), and the climatic water deficit (from TerraClimate, Section 2.5) during the active season ranged from below 10 to above 50 mm (Fig. 1). The forests are evergreen needleleaf (ENF, 7 forests), deciduous broadleaf (DBF, 6), and mixed (1) which we analyzed with the DBF group. Most sites are natural vegetation except CA-TP3 and CA-TP4 which were originally planted as monocultures (in 1974(in and 1939(in , respectively, Arain et al., 2022)).
Flux tower data were primarily obtained from the AmeriFlux database (https://ameriflux.lbl.gov/). The R package REddyProc (version 1.3.2) was used to remove periods of low turbulence using a site-specific friction velocity threshold, to gap-fill NEE and weather data (Wutzler et al., 2018), and to partition NEE (REddyProc variable NEE_U50_f) into GPP (GPP_U50_f) and ecosystem respiration using the nighttime method (Reichstein et al., 2005). We avoided the daytime method of Lasslop et al. (2010) due to erroneous GPP in winter (Bowling et al., 2024), which would lead to inaccurate determination of the timing of seasonal transitions. Failure of the algorithms to identify suitable friction velocity thresholds occurred in a few cases, leading to entire site-years failing the REddyProc gap-filling process. Switching to the FLUXNET2015 database (https://fluxnet.org/data/fluxnet2015-dataset/) for US-UMB and the AmeriFlux FLUXNET product (https://ameriflux.lbl.gov/data/flux-data -products/oneflux-processing/) for CA-Cbo and CA-TP3 alleviated this site-dependent problem and enabled the use of longer records at these sites.
Half-hourly GPP was summed to calculate total GPP during particular seasons (Section 2.3). Because cumulative GPP is sensitive to data gaps, years with unfilled gaps during the active season were not included, and data from all REddyProc quality flag categories were used. Years with high cumulative GPP quality flag sums (>1000) indicating poor data quality were removed. In addition, we found no association between the quality flag sum and outliers of cumulative GPP during the active season (see list of Abbreviations). Processed 30-min flux data were de-spiked by binning half hourly data into 13-day windows and identifying data above or below the median ± 4x the median absolute deviation, separately for day and night (Papale et al., 2006). These 1 and Appendix B for more site details.
spikes were replaced by the average GPP of that day before calculating seasonal sums.
The method of Bowling et al. (2024) was used to determine the timing of seasonal transitions for photosynthesis. Full details and code to calculate transitions may be found in that paper. Briefly, the method evaluates the response of NEE to photosynthetically active radiation (PAR) in 5-day moving windows to calculate the seasonal pattern of whole-forest photosynthetic capacity at high light, which we refer to as GPP 1800 . This quantity is the value of a fitted curve (not shown) between NEE and PAR during each 5-d window at a PAR level of 1800 µmol m -2 s -1 , after adjusting for respiration. The annual pattern of GPP 1800 is shown in Fig. 2a for one site (US-Ha1). The annual GPP 1800 time series were fitted with 2 logistic equations (not shown), and the 10 and 90% thresholds between baseline and summer maximum of the logistic fits were used to define dates of transition between seasons (SOS and EOS at 10% threshold, and transitions with the active season at 90%, Fig. 2a). We define the active season (AS) for GPP as the time period between SOS and EOS (this is the main period of carbon uptake), and further divide this into three periods (spring ramp, summer, and autumn ramp) based on the 90% threshold crossings.
To test hypotheses, we calculated cumulative GPP (g C m -2 ) in each portion of the active season (full active season, spring ramp, summer, autumn ramp). Years with missing seasonal transition dates (SOS, EOS, etc.) due to missing data or poor-quality logistic fits prevented seasonal identification and were excluded from the analysis. The active season length was defined as the number of days between the SOS and EOS, inclusive. Cumulative GPP in each season is referred to as ΣGPP with a subscript indicating season (ΣGPP AS , ΣGPP spring ramp , etc.).
The Snow Data Assimilation System (SNODAS) is a data-constrained reanalysis product that combines satellite, airborne, and ground data with models of weather prediction and snow energy and mass balance (Barrett, 2003). Daily, 1km 2 gridded SNODAS data were obtained (accession date June 15, 2022) from the National Snow and Ice Data Center (National Operational Hydrologic Remote Sensing Center, 2004). For sites in the contiguous US, SNODAS data were available from 2004present; for sites in southeast Canada, data were available from 2010present. To reduce the impact of random uncertainty associated by using a single pixel for each flux tower, we averaged SNODAS precipitation (solid and liquid) and SWE variables for all pixels contained or partially-contained within a 2 km radius of the flux tower with similar Fig. 2. Overview of our phenological framework and associated seasonal definitions, using data for the 2013 active season and prior dormant season from Harvard Forest (US-Ha1). a) Time series of the light response of photosynthesis, evaluated at high light (GPP 1800 ), were analyzed, using the method of Bowling et al. (2024), defining start of season (SOS, pink circle), start and end of summer (yellow circles, binding the yellow box highlighting "summer"), and end of season (EOS, pink circle). The period between SOS and the start of summer defines the "spring ramp" (green box). The "autumn ramp" (gold box) is the period between the end of summer (yellow circle) and EOS (pink circle). The "active season" for photosynthesis is the combination of spring ramp, summer, and autumn ramp periods (top olive box). The "dormant season" is the cold-season period between EOS in autumn and SOS in the subsequent spring (blue box). b) Time series of snow water equivalent (SWE), directly observed at US-Ha1 (gray) and from the SNODAS SWE model product (blue). The DoY and magnitude of peak SWE and the date when snow has fully disappeared (snow all gone, SAG) were obtained from the SNODAS data (arrows).
vegetation cover class, based on MODIS IGBP land cover. We overlaid these pixels with a digital elevation model (Amante and Eakins, 2009) before taking the average weighted by similarity in elevation to the flux tower, where each pixel's weight = 1/ absolute value of (elevation of pixelelevation of tower). Alternate buffer sizes and weighting schemes, as well as filtering for outliers based on slope, aspect and leaf area index were also considered. We selected the final processing method after comparing how these variations influenced comparison with in-situ SWE data at three sites where snowpack data were available: US-Ha1, US-NR1, and US-GLE (see Appendix A). The latter site was used for SNODAS evaluation but not included in GPP analyses due to major insect disturbance (Frank et al., 2014).
The SNODAS data were used to calculate metrics of the magnitude and timing of snowmelt for each flux tower. Date of snow disappearance, or snow all gone (SAG), was defined as the first day after which no new SWE was present (see Fig. 2b for a comparison of SNODAS and in-situ SWE hydrographs). Comparison of SAG determined from SNODAS versus in-situ data show a slight overestimation of SNODAS SAG at US-NR1 and US-GLE, and at US-Ha1 a few years exhibited discrepancy between SNODAS and in-situ SAG (Fig. A2a-c). This discrepancy was caused by years where snowpack disappearance was followed by a few small, isolated snowfall events not present in the SNODAS product (Fig. A2d), however the SNODAS product appears to give a good indication of when the primary snowpack has disappeared. Total snowfall (mm) was calculated as the sum of the SNODAS variable solid precipitation (PRSL). Total rain (mm) was calculated for each season as the sum of liquid precipitation (PRLQ). Peak SWE (mm) was determined as the annual maximum SWE from SNODAS (SWEM variable), and the timing of peak SWE as the day of year (DoY) on which it occurred. The length of the melt period was defined as the number of days between peak SWE and SAG, inclusive. Solid fraction (percentage of annual precipitation falling as snow) was calculated as 100xPRSL/(PRSL+PRLQ).
Soil moisture data were used to examine water availability for plants. All available soil volumetric water content (VWC) data for each forest were accessed from AmeriFlux, or obtained directly from site scientists. For sites that included observations across a soil depth profile, we assessed whether the use of profile-integrated VWC significantly affected our results compared to the use of single depth measurements. We found that there was no improvement with integrated profile measurements, and therefore in favor of consistency across sites, we used VWC data from a depth of 15 cm which were available at all but 1 site. VWC data were not available for all years at all sites, and were not available for US-Vcm.
The climatic water deficit (CWD) was used to examine the combined effects of water availability and vapor pressure saturation deficit of air on moisture limitation for plants, obtained from TerraClimate (Abatzoglou et al., 2018) at 4 km resolution. The CWD is calculated as the difference (mm) between reference evapotranspiration and actual evapotranspiration. Reference evapotranspiration is calculated by Ter-raClimate assuming standard parameters for a grass surface (e.g., Allen et al., 1998), which can be problematic applied to forests and is likely overestimated (Sun et al., 2016). Actual evapotranspiration is calculated by TerraClimate using a Thornthwaite-Mather water balance model (Dobrowski et al., 2013), with additional uncertainty. However, the CWD is biologically meaningful (Stephenson, 1998) and has been shown to be a robust metric of plant water relations in studies of productivity and forest mortality (e.g., Anderegg et al., 2015;Hoylman et al., 2019). The CWD data were available monthly, so to match the timing of seasonal transitions, we calculated the average CWD of all months contained within a season (active season, spring ramp, etc.), weighting each month by its proportion contained in the season (e.g. if SOS occurred on April 25, then April CWD was weighted to be 5/30ths of the overall mean spring ramp CWD).
Linear regression, correlation analysis, and mixed effects models were used to test hypotheses. To test H1, we examined simple linear regressions between active season length and SOS, SAG, and EOS. To test H2, we regressed ΣGPP AS versus active season length and SAG. Correlation analysis was used to quantify bivariate relations between cumulative GPP (ΣGPP AS ) and active season length, active season CWD, mean T air during active season, length of snowmelt, peak SWE, SAG, total snowfall, spring rain, and summer rain (to test H2). We also examined correlations between cumulative GPP during the spring and autumn ramp and summer seasons with cumulative rainfall in each season (providing moisture-based alternatives that might help explain H2).
Mixed effects models were further used to test hypotheses, constructed using a top-down model selection process (Zuur et al., 2009), which uses iterative backwards selection to find those models that explain the most variation with the minimum necessary parameters. The ENF and DBF forest types were analyzed separately. Separate models were built with ΣGPP AS and active season length as response variables, and applied separately for each hypothesis.
To test H1, we built linear mixed effect models separately for ENF and DBF with active season length as the response variable. First, we fit a saturated fixed-effects-only model with all possible terms that represent biologically real hypotheses, and their interactions. For H1, starting variables included: mean T air , PAR, and VPD for each season (spring ramp, summer, autumn ramp), total rainfall in each of these and active season, timing of SAG, amount of peak SWE, timing of peak SWE, total snowfall, mean VWC for each season, timing of spring maximum (VWC max ), summer VWC min , and mean active season CWD. We used the dredge function in the MuMIn package in R (Bartoń, 2023) to determine the relative importance of each candidate fixed effect variable based on the ranked Akaike information criteria (AIC). This estimate of variable importance is made by summing the AIC weights across all candidate models which contain that variable (Burnham and Anderson, 2004). Fixed effects variables with high importance (> 0.8) were then used to construct a less saturated model with both 'site' and 'year' considered as possible random effects. Next, we used the step function in the lmerTest package in R (Kuznetsova et al., 2017) to perform backward elimination of random-effect terms followed by backward elimination of fixed-effect terms to find the most parsimonious model, as follows. First, random effects (site and/or year) were eliminated based on the likelihood ratio test. In all cases, the best random effect structure was a random intercept model with 'site' as a random effect. Random slopes led to overfitting in all cases and were therefore not included. Then, fixed effects were eliminated based on ANOVA with p-values calculated using Satterthwaite's method (Kuznetsova et al., 2017). In some cases, the selected model resulted in singular gradient errors or convergence failures indicating overfitting; for these cases, fixed effect variables were dropped one at a time using likelihood ratio tests of nested models to determine the final model. All final models were tested for collinearity among independent variables, such that the variance inflation factor was < 2 for all retained variables. All model comparisons were made using maximum likelihood (ML) fitting, then final models were presented using restricted maximum likelihood (REML) estimation, with marginal and conditional R 2 calculated using Nakagawa and Schielzeth's (2013) method for mixed models.
To test H2, we followed the model selection method above to determine the most important explanatory variables for ΣGPP AS and build parsimonious models for ENF and DBF. For H2, starting candidate variables included: the timing and magnitude of peak SWE, total snowfall, SAG, total rainfall in each season, mean T air , PAR, and VPD in each season, length of spring ramp, SOS, active season CWD, active season length, mean active season VWC, and the days of year of spring VWC max and summer VWC min .
In addition, we tested H2 and H3 spatially using mean annual conditions for each site rather than values subject to interannual variability. All data from both ENF and DBF were combined to construct the model (including biome as a fixed effect did not provide any significant improvement based on likelihood ratio testing of AIC). We tested H2 using a mixed effects model with ΣGPP AS as the response, and the threeway interaction between active season length, mean annual summer CWD, and mean annual solid fraction. Other interannually-averaged candidate variables that were considered but did not show improvement (based on likelihood ratio testing of AIC) were MAT, MAP, mean active season precipitation, mean annual peak SWE, and mean active season CWD. We then used the same approach for H3, where SAG was used in place of active season length-with total ΣGPP AS the response, and SAG, mean annual summer CWD, and mean annual solid fraction as fixed effects.
Interannual variation in active season length was significantly correlated with SOS (the date when GPP first reached 10% of maximum photosynthetic capacity) at most of our study forests (Fig. 3a, Table 2). This is not surprising as SOS and EOS define the active season length. This pattern was also present when all sites were analyzed together (a single regression combining all sites in Fig. 3a was highly significant with R 2 of 0.87). However, the active season length was not significantly correlated with EOS at any individual sites, though they were correlated with sites combined (Table 2). These results indicate that the initiation of photosynthesis in spring was the primary determinant of variability in active season length. This is a necessary requirement supporting H1 (variation in active season length is related to variation in SOS), but full support for H1 requires linkage between active season length and date of full snow disappearance (SAG). Regressions of active season length with SAG from SNODAS generally had negative slope (Fig. 3b, testing H1), but the correlations were weak and slopes were not significantly different from zero (Table 2, CA-TPD was an exception). Regressions of SOS with SAG were also weak and mostly non-significant (Table 2), except when all sites within a forest type were analyzed together. These results do not support H1.
Bivariate correlation analysis was used to examine how active season length and environmental variables influenced interannual cumulative GPP at each site (Fig. 4). Many significant correlations were present, both positive and negative, and were highly site-specific. In general, rainfall in each season was correlated with the GPP for that season, but clear patterns were not present at all sites. There were significant negative correlations between timing of snowmelt (SAG) and ΣGPP AS at 2 sites (US-Ha1, US-Ha2), indicating earlier melt led to higher productivity (some support for H3). In other sites, that relationship was not significant (not supporting H3). The ΣGPP AS was significantly correlated with active season length at 3-4 sites (US-Ha1, US-Ho2, and US-Vcm at p < 0.05, US-Ho1 at p < 0.1) but not others (mixed support for H2). Peak SWE positively influenced ΣGPP AS at US-Wcr only, correlations for other sites were not significant (general lack of support for H3). There did not appear to be uniform consistency in the direction or strength of relationships due to biome or site average active season CWD. There was strong correlation between ΣGPP autumn and autumn rain (Fig. 4), perhaps due to warmer autumns being wetter (anomalies of T air and GPP were both positively and significantly correlated with autumn rain, data not shown).
Overall, the data in Fig. 4 indicate that the environmental drivers of active season cumulative GPP differed by site, and GPP was in general not well characterized by interannual variation in weather. Therefore, we performed cross-site analyses to test the hypotheses, and present them separately here. We employed a model selection approach using linear mixed effects models. The inclusion of site as a random effect 2. accounted for the large degree of inter-site variation that was not well characterized by phenological or climatic drivers, such as species composition, forest age, soil type, nutrient limitations, etc.
To assess the importance of snowmelt and other climatic drivers on active season length, linear mixed effect models were built separately for ENF and DBF with active season length as the response variable. Starting candidates included: mean T air , PAR, and VPD for each season (spring, summer, fall), total precipitation in each season and the entire active season, timing of SAG, amount of peak SWE, timing of peak SWE, total snowfall, mean VWC for each season, timing of spring maximum (VWC max ), summer VWC min , and mean active season CWD.
Final models selected with standardized coefficients are shown in Table 3, and their graphical representation is shown as marginal effects plots in Fig. 5. Marginal effects plots show the partial residuals for each fixed effect term after holding all the other terms constant at their median. One can consider the X and Y axes of a marginal effects plot as 'X and Y after all other predictors from the model have been accounted for,' and the slope of each line in Fig. 5. represents the partial regression coefficients (see Table 3). In addition, in Table 3 it is useful to compare the variance explained by fixed effects (marginal R 2 ) with the variance explained by both fixed and random effects (conditional R 2 ).
For the ENF biome, variables with high importance (> 0.8) were the date of peak SWE, SAG, mean spring T air , mean autumn PAR, total autumn rain, and AS CWD. (See Appendix B for information about the interannual variability of parameters important in the mixed model results.) The selected model included (in order of standardized effect size) mean autumn PAR, total autumn rain, the timing of SAG, and the timing of peak SWE (Table 3, marginal R 2 =0.75, conditional R 2 =0.85). That AS length declined with timing of SAG and timing of peak SWE provides some support for H1 for ENF. It appears that, in addition to the timing of snowmelt, autumn conditions have a degree of control over AS length after partial pooling of sites, though this was not true at individual sites (see Fig. 4 and AS length vs EOS in Table 2).
For the DBF biome, variables with high importance were mean spring and autumn T air , total spring and autumn rain, mean VWC autumn , and timing of spring VWC max . The final model was built, in order of effect size, using mean spring T air , total autumn rain, and timing of spring VWC max (Fig. 5, Tables 3, B1). Longer active seasons had cooler springs on average (Fig. 5a). Similar to ENF, longer active seasons were associated with more autumn rain (Fig. 5c,e). The marginal explanatory power was low and the conditional power was high (R 2 = 0.32 and R 2 = 0.85, respectively), indicating that the random effect due to site accounted for more variance than the fixed effects. Thus, the length of the active season in DBF was poorly constrained by these environmental predictors, despite being the most parsimonious model found. This does not support H1 for DBF. That warmer springs were associated with shorter active seasons does not necessarily say that warming shortened the active season, but simply reflects that later SOS results in warmer spring temperature. Rather, it would be expected that warmer air will increase the length of the active season in DBF (Baldocchi et al., 2018). As a caveat, note that the length of seasons varies between years and sites, and thus variables that are summed (such as autumn rain here) are somewhat problematic to interpret as a result.
We followed the above approach for H2 to determine the most important explanatory variables for ΣGPP AS (as response variable), separately for ENF and DBF. Starting candidate variables included: active season length, the timing and magnitude of peak SWE, total snowfall, SAG, total rainfall in each season, mean T air in each season, length of spring ramp, SOS, active season CWD, mean active season VWC, DoY of spring VWC max , and summer VWC min . For the ENF biome, variables with high importance included active season length, SAG, mean active season T air , and active season CWD. The final model was built using active season length and active season CWD (Fig. 6b,c, Tables 4, B1). For the DBF biome, variables with high importance were active season length, the timing of peak SWE, mean active season VWC, mean active season T air , and the length of spring ramp, however the final model included only the effect of active season length (Fig. 6a, Tables 4, B1).
For both ENF and DBF, active season productivity increased with the length of the active season (Fig. 6a,b), providing support for H2 for both forest types. The effect of active season length was larger in DBF compared to ENF (3.9 ± 1.8 gC m -2 d -1 versus 2.0 +-0.8 gC m -2 d -1 , respectively), and in ENF ΣGPP AS also declined with increasing active season CWD (Fig. 6c, Table 4). For comparison, Launiainen et al. (2022) found increasing trends in GPP with lengthening active season over many years (~ 8 g C m -2 year -1 ), and Baldocchi et al. (2001) reported a general interannual pattern of higher GPP with longer active season of 5.7 g C m -2 d -1 across temperate DBF sites, which were not constrained to seasonally-snow covered sites. For both ENF and DBF models, however, the explanatory power of these fixed effects was very low (marginal R 2 =0.11 and R 2 =0.06, respectively), while site as a random effect explained most of the variance (conditional R 2 = 0.90 and R 2 =0.72, respectively; Table 4).
Mixed effects model predictor coefficients for the response variable active season length (days) in both ENF and DBF (see Fig. 5). The coefficients are standardized to compare the magnitude of fixed effects within the model. Values in parentheses are standard errors. τ 00 is the random effects variance, and σ 2 is the model residual variance. Marginal R 2 is the variance explained by fixed effects, while conditional R 2 is the variance explained by both fixed and random effects. p-values * (p < 0.05), ** (p < 0.01), *** (p < 0.001).
ENF DBF mean PAR autumn -0.132*** (0.033) total autumn rain 0.066*** (0.016) 0.049*** (0.014) DoY SAG -0.293** (0.104) DoY peak SWE -0.153** (0.052) mean spring T air -4.438*** (0.675) DoY spring VWC max -0.246** (0.083) Groups: Site 7 7 τ 00 103.8 214.48 σ 2 152.4 63.24 marginal R 2 0.75 0.32 conditional R 2 0.85 0.85 Observations 77 66 Akaike Information Criterion 642.626 496.864
Next, we evaluated whether variation in seasonal precipitation influenced the linkage between active season length and ΣGPP AS . Through the above correlation analysis between active season length and GPP within individual sites (Fig. 4), as well as model selection within biomes with partial pooling across sites, we were unable to find consistent patterns that supported or refuted H2, due to the large amount of unexplained interannual variability within and across sites (Figs. 5,6). Therefore, we tested H2 spatially, using mean annual conditions for each site rather than including interannual variation.
Shown in Fig. 7 (and Table 5) are linear regressions of ΣGPP AS with active season length and SAG, as a function of mean summer CWD averaged across all years for each site (color axis). Sites that had a higher mean annual summer CWD tended to have a less positive, or even negative, GPP response to active season length (Fig. 7a) and a less negative, or positive, response to SAG (Fig. 7b). These results indicate that summer moisture deficit potentially explains some of the variation in the linkage between ΣGPP AS and active season length (caveats to allow support for H2). Since water from the winter snowpack might influence summer soil moisture availability for plants (via the moisture effect), and thus influence summer CWD, we built a mixed effects model with ΣGPP AS as the response variable, and the three-way interaction between active season length, mean annual summer CWD, and mean annual solid fraction. Other interannually-averaged candidate variables that were considered but did not show improvement were MAT, MAP, mean annual active season rainfall, mean annual peak SWE, and mean annual active season CWD. Data from both ENF and DBF biomes were combined to construct the model, and including biome as a fixed effect did not provide any significant improvement based on likelihood ratio testing of AIC. Partial effects of the three-way interaction are shown in Fig. 8. All terms and interactions in the model were statistically significant (p < 0.05) and it was the most parsimonious combination variables (lowest AIC) based on nested likelihood ratio tests.
The model (Fig. 8) significantly improved the prediction of ΣGPP AS (marginal R 2 = 0.50, conditional R 2 =0.88) compared to the previous models (Fig. 5, Table 3) that did not account for site average climate conditions. For sites with higher mean annual CWD (those with drier summers on average), GPP decreased with longer active seasons (red lines in Fig. 3.8). For sites with wetter summers on average, GPP increased with longer active seasons (blue lines in Fig. 3.8). The magnitude of these effects decreased as solid fraction increased, and ΣGPP AS became less dependent (shallow slopes) on the active season length. These results provide mixed support for H2, based on CWD.
To test H3, we used a similar approach as for H2, with SAG used in place of active season length. We used ΣGPP AS as the response variable, and SAG, mean annual summer CWD, and mean annual solid fraction as candidate fixed effects. Rather than a three-way interaction, however, only the interaction between SAG and summer CWD was significant (p <
Mixed effects model predictor coefficients for the response variable ΣGPP AS (g C m -2 ) in both ENF and DBF (see Fig. 6). Coefficients are unstandardized to facilitate comparison of the fixed effect size of active season length between the two models. Values in parentheses are standard errors. τ 00 is the random effects variance, and σ 2 is the model residual variance. Marginal R 2 is the variance explained by fixed effects, while conditional R 2 is the variance explained by both fixed and random effects. p-values * (p < 0.05), ** (p < 0.01), *** (p < 0.001).
ENF DBF active season length 2.027** (0.751) 3.921* (1.808) active season CWD -3.214** (1.245) Groups: Site 7 7 τ 00 86,887 74,666 σ 2 10,518 31,838 marginal R 2 0.11 0.06 conditional R 2 0.90 0.72 Observations 77 69 Akaike Information Criterion 954.750 925.053 0.05). The partial residuals of the interaction between SAG and summer CWD are shown in Fig. 9, with the partial residuals of mean annual solid fraction shown in the inset. The explanatory power of the model was similar but slightly lower than the corresponding model with active season length (marginal R 2 =0.45 and conditional R 2 =0.86). These results demonstrate that in areas with low summer CWD (favorable moisture conditions), later snow disappearance has a negative effect on active season GPP, but this relationship diminishes with increasing CWD. Peak SWE never exhibited high importance in our models. These results provide mixed support for H3, based on CWD.
Studying ecosystem dynamics on interannual scales is challenging conceptually and statistically. Conceptually, there are many environmental and biotic drivers of interannual variability in C fluxes operating at different timescales with unknown lags and legacy effects (Richardson et al., 2010). However, it is difficult to distinguish between the independent effects of particular drivers on carbon fluxes. The sensitivity of carbon fluxes to climatic variability appears to progressively decline or become more difficult to detect at increasing timescales, and our results support the general consensus in the literature that interannual variation in GPP is not well-represented by direct responses to fluctuating environmental conditions (De Pue et al., 2023;Richardson et al., 2007;Stoy et al., 2009;Urbanski et al., 2007;J. Wu et al., 2012). GPP and other ecosystem processes respond to environmental variation in changing ways throughout the annual cycle (Launiainen et al., 2022). Statistically, it is difficult to obtain records of sufficient length to find significant relationships, as interannual variability in ecosystem C fluxes is large relative to any trend driven by a specific environmental variable (Baldocchi et al., 2018). Although the lengths of many eddy covariance records are now multi-decadal, the statistical significance of reported relationships between anomalies in active season length and carbon fluxes often remains weak (Richardson et al., 2009).
We tested three hypotheses, with mixed results for all. We found that a longer active season for GPP occurred when photosynthesis started earlier in the year, but interannual variation in GPP was not directly affected by timing of snowmelt at individual sites (Fig. 3, Table 2). When sites were considered together with the use of mixed effects models, results indicated that later snowmelt timing had a negative effect on total active season length for ENF, but was not important for DBF (Fig. 5f, Table 3). We found significant correlations between environmental variables and ΣGPP AS at some sites, but not others, and there was large variability in direct predictors of interannual GPP even among sites that were in the same functional class and climate zone, and relatively close geographically (Fig. 4). Mixed effects models indicated that, when sites were combined, ΣGPP AS was significantly and positively affected by active season length (Fig. 6) in both forest types. Further exploration with mixed effects models indicates that this pattern is dependent on long-term average CWD and solid fraction (Figs. 7,8,9). These results, taken together, highlight the complicated nature of the controls on interannual variability of forest carbon sinks. Our attempt to find direct predictors of interannual variation of GPP within sites was limited. Relationships that were found differed from site to site with no evident patterns between sites. That relationships were highly site-specific is consistent with other studies (see review by Baldocchi et al., 2018). Evidence suggests that recent warming trends have led to earlier SOS over the last few decades (Badeck et al., 2004;Jiang et al., 2023;Richardson et al., 2006). It is well established that delayed spring onset results in shorter active seasons for GPP (C. Wu et al., 2012b), and in the present study this was true at almost all sites (Fig. 3), capturing both interannual variability, as well as spatial (R 2 = 0.87). We did not however find that autumn senescence was correlated with interannual variability in the active season length at individual sites (Table 2), in contrast with previous studies (Desai et al., 2022;Fu et al., 2017;Keenan and Richardson, 2015;C. Wu et al., 2012a). At US-NR1, interannual variation in the active season length was determined more by the duration of snow melt than the timing of senescence (Monson et al., 2005).
We did not find that peak SWE or timing of full melt (SAG) were useful to explain interannual variability in active season length using simple correlation at each site. However, SAG and timing of peak SWE were the most important variables for explaining interannual variability in the length of the active season across sites in ENF (Fig. 5). This was not true for DBF, however the inclusion of spring VWC max timing suggests some importance of spring hydrology (Fig. 5). Overall, active season length in DBF was less responsive to environmental drivers with a low marginal R 2 (see Table 3) indicating that fixed effects explained a very low proportion of variability. This perhaps reflects the general understanding that while ENF active season length is contingent on the seasonal relief of environmental constraints (Bowling et al., 2018;Monson et al., 2005), in DBF the carbon uptake period can also be constrained by new leaf emergence (Barr et al., 2007;Desai et al., 2022;Gu et al., 2003).
The active season length was the most important predictor for variation in ΣGPP AS for both ENF and DBF across sites, despite only being a significant predictor of interannual variation at three sites (Fig. 4). Similarly, Wu et al. (2012b) found that active season length was a good indicator of spatial variability in annual net ecosystem productivity of North American forests, but that predictors which had strong spatial correlation were not good predictors of interannual variability. Numerous past studies have demonstrated positive relationships between NEE and active season length in ENF (Danielewska et al., 2015), DBF (Baldocchi et al., 2001;Desai et al., 2022;Finzi et al., 2020;Gu et al., 2003;Richardson et al., 2010Richardson et al., , 2009;;White et al., 1999), or both (Churkina et al., 2005;Fu et al., 2017). In agreement, our results showed a positive association in both biomes (Fig. 6), and the effect of prolonged active season length on ΣGPP AS was stronger in DBF than ENF (Table 4). Our model selection suggested that this could be attributed to the mediating negative relationship between GPP and CWD in ENF (Fig. 6b, c). Because CWD accounts for both reference and actual evapotranspiration, it has the benefit of integrating over both the effects of soil moisture supply and atmospheric water vapor demand, which are often correlated and difficult to disentangle. However, the CWD is not a measure of actual stand water use with respect to water availability, and ignores species and stand-level controls on evapotranspiration (e.g., Fu et al., 2022;Launiainen et al., 2016). That CWD was better suited in the role of mediator compared to snow, rain, or soil moisture metrics alone, agrees with previous studies that have highlighted the importance of the vapor pressure deficit to limiting canopy conductance in mesic forests (Novick et al., 2016). Regardless, we found the relationship between active season length and ΣGPP AS had very weak marginal explanatory power for both ENF and DBF, while the site-level variation (explained by random intercepts in the statistical model) was extremely dominant (compare marginal and conditional R 2 , Table 4). These results highlight the large degree of variability between sites that is not well represented by seasonal or annual weather.
In all instances where snowpack variables were important in determining the active season length and ΣGPP AS , there was a greater importance of snow timing metrics compared to the magnitude of peak SWE or total snowfall amount, in agreement with Knowles et al. (2018). This included 1) the importance of SAG in determining the active season length in ENF (Fig. 5c); 2) the appearance of timing metrics as statistically important for ΣGPP AS in ENF and DBF (though they did not end up in the final models); and 3) the finding that SAG was related to ΣGPP AS across sites after accounting for site differences in mean annual precipitation solid fraction and summer CWD (Fig. 9). We would expect that snowmelt timing would be an important factor for active season length and ΣGPP AS , under the premise that in some locations soil moisture is highest in spring and declines after snow is fully gone. In western U.S. conifer forests it was shown that peak annual soil moisture coincides with the date of snow disappearance (Harpold et al., 2015;Harpold and Molotch, 2015), yet we found this to be true only at the snow-dominated high-elevation subalpine forest site US-NR1 (R 2 =0.69, p < 0.001, data not shown). And although all sites showed a decline in mean VWC from spring to summer, mean spring VWC did not show significant increase compared to mean winter VWC at 8 of the 13 sites with available soil moisture data (not shown). These sites did not follow the textbook hydrologic dynamics of seasonally snow-covered forests that have stable dormant season soil moisture followed by a clear spring melt period (e. g., Maurer and Bowling, 2014), but rather experienced influxes of snowmelt and rain during the winter. These patterns complicate attempts to understand the impact of snowmelt timing at these sites, and likely contribute to the variability in the interannual explanatory power of SAG at individual sites (Fig. 3, A2). Nevertheless, it appears that the timing of peak SWE and SAG is still important when considered across sites, despite the lack of coincidence with peak soil moisture. We did not find that peak SWE was important for GPP. At the site level, only US-WCr had a significant relationship between interannual ΣGPP AS and peak SWE (Fig. 4), as previously documented by Desai et al. (2022), though they attributed this to a soil temperature rather than moisture effect. Wang et al. (2018) suggest that the dependence of summer productivity on snowmelt is determined by both the legacy effect of winter SWE on active season soil moisture and by the degree to which vegetation growth is water limited. Regarding the former, there were mixed, inconsistent results as to whether there was a legacy effect of the amount of SWE on summer soil moisture at the site level (data not shown). For instance, total solid precipitation was significantly correlated with minimum summer VWC at US-NR1 and CA-Cbo, however we did not find that this translated into a legacy effect of total snow on active season or summer GPP at these sites (Fig. 3). Similarly, Richardson et al. (2009) found that the lagged effect of spring phenological anomalies on summer fluxes was weak and non-significant, due to the larger influence of summer weather. Alternatively, there may be important seasonal moisture legacies associated with land-atmosphere teleconnections that we are missing. For example, high snow years in the Rocky Mountains are associated with lower North American Monsoon rainfall in the subsequent summer (Lo and Clark, 2002;Notaro and Zarrin, 2011). Patterns of seasonal water use by trees differ across the western US due to the spatial gradient in monsoon rainfall (Szejner et al., 2016).
A primary goal was to determine whether we could detect which environmental controls determine the relative dominance of the moisture vs. growth period effects of earlier snowmelt. We found that after accounting for the interannual variability at each site, mean annual summertime CWD and mean annual solid fraction mediated the response of ΣGPP AS to variation in active season length and SAG, and this was not dependent on forest type (Fig. 8). For sites with higher mean summer CWD (those with drier summers on average), ΣGPP AS decreased with longer active seasons (red lines in Fig. 8). For sites with wetter summers on average, ΣGPP AS increased with longer active seasons (blue lines in Fig. 8). The magnitude of these effects decreased as solid fraction increased, and active season GPP became less dependent on the active season length (Fig. 9). In forests with low summer CWD (sufficient moisture), later SAG had a negative effect on active season production, and this relationship diminished with increasing CWD (Fig. 9). These average climatic conditions were found to be more important in determining the magnitude and direction of the relationship between active season length and ΣGPP AS than when CWD or the annual precipitation solid fraction were examined interannually (Fig. 3). Thus, the opportunity for forests to capitalize on the C production potential of longer active seasons appears to be dependent on the degree of reliance of vegetation to average moisture conditions determined by long-term climate characteristics rather than interannual fluctuations in weather. For natural vegetation this may be a result of competition and adaptation to local microclimate. This highlights the importance of considering longer-term ecological and demographic processes when trying to predict how vegetation will respond to future changes in climate. Previous studies have also found that average factors acting over long time scales, such as water table depth (Dunn et al., 2007), mean annual temperature (White et al., 1999), or the average vertical distribution of soil moisture (Martin et al., 2018), mediate sensitivity to changes in active season length. Our results add to the growing body of literature which has shown that drier sites are vulnerable to increasing summer drying in response to longer active seasons (Knowles et al., 2018;Parida and Buermann, 2014), and that this is true across different biomes (Buermann et al., 2018;Butterfield et al., 2020;Xu et al., 2020).
Some limitations apply to our study. First, our hypotheses are almost certainly too simplistic, particularly given the seasonal nature of rain and snowfall. Patterns of forest productivity in response to seasonal hydrologic variation, and its future change, are likely to vary for different seasonally snow-covered climates, such as those influenced by large-scale continental patterns of precipitation in the western US (Trujillo and Molotch, 2014), Europe (Beniston et al., 2018), and the Asian and North American monsoons (Adams and Comrie, 1997;Wu and Qian, 2003). Second, present availability of snowpack and related moisture data at flux towers is quite limited. We recommend that flux tower scientists in seasonally-snow-covered biomes consider the addition of continuous electronic snow depth and SWE instrumentation, and snow temperature (which indicates the melt process, Burns et al., 2014) as a part of the standard suite of environmental observations. We also share the building enthusiasm to include observations of soil and plant water potential, which are among the most useful metrics of plant physiological response to water availability (Novick et al., 2022). Third, geographical representation of SNODAS products severely limited the number of flux towers used and their representation across climate space (most of our study forests were similar in climate, Fig. 1, and they were all in North America). Fourth, the spatial resolution of the gridded SNODAS (1km 2 ) and TerraClimate (4km 2 ) products is coarse relative to flux tower footprints (Chu et al., 2021). These data-availability limitations were a constraint on the breadth to which we were able to test our general hypotheses. Finally, we acknowledge that we have ignored uncertainty in the standardized GPP products. In general there is large interannual variation in NEE, GPP, and ecosystem respiration, when compared to the multi-year mean at any site (Baldocchi et al., 2018), and there is likely to additional variation among the many GPP products available (Pastorello et al., 2020). Future analyses that examine the growth period and moisture effects in the context of climate and environmental change will be strengthened if we can alleviate these challenges.
We synthesized 145 site years of eddy covariance flux data from 14 deciduous and evergreen forest sites in the US and southeast Canada with gridded SWE and precipitation estimates from SNODAS. We used these data to study the spatiotemporal response of active season GPP to interannual and spatial variability in active season length, timing of snowmelt, and the date of disappearance of snow. We found that the relative dominance of the moisture and growth period effects of earlier snowmelt and associated longer active seasons was dependent on siteaverage moisture conditions, determined by long-term summer climatic water deficit and precipitation solid fraction, rather than interannual fluctuations in weather. In addition, we did not find that the magnitude of peak SWE was important for determining the active season length or total GPP. SWE did not appear to have a legacy effect that influences active season GPP. However, lagged effects may have been overshadowed by the more direct influence of active season weather. Finally, we emphasize that although our hypotheses were supported across broader spatial and temporal scales, they were not correspondingly supported within individual sites. There was a large degree of interannual variability both within and between sites that was not well represented by seasonal climate or phenology. anomalous data points at US-Ha1 were found and illustrated in detail in d) which shows that a year with snowpack disappearance followed by a few small, isolated snowfall events can in a discrepancy in the SAG estimate.
Interannual variability of active season cumulative GPP and active season length, important dates, and variables found to be important for the mixed models in Tables 3 and 4. Data shown as the interannual mean and standard deviation (in parentheses). Actual years analyzed were smaller than the full range available for some sites, based on data availability (particularly SNODAS and flux data overlap) and quality assurance. Exact years analyzed can be found in the supplemental dataset. Mean annual peak SWE and mean annual solid fraction were calculated from SNODAS, all other data were obtained from the literature cited. Sites are ordered by mean annual temperature for each forest type as in Table 1.
Site Biome Number of Years Analyzed ΣGPP AS (g C m - 2 ) AS length (d) SAG (DoY) SOS (DoY) EOS (DoY) peak SWE (mm) DoY of peak SWE (DoY) Mean Annual Solid Fraction (%) Mean Spring T air ( • C) DoY of Spring Max. VWC (DoY) AS Climatic Water Deficit (mm) Mean Autumn PAR (μmol m -2 s -1 ) Sum Autumn Rain (mm) US-Syv DBF 8 898 (127) 178 (9) 123 (9) 124 (9) 248 (14) 172 (58) 69 (20) 30 (0.05) 12.9 (1.9) 147 (9) 6.9 (4.5) 252 (32) 130 (84) US-WCr 10 883 (122) 144 (8) 120 (13) 135 (5) 245 (10) 116 (49) 63 (21) 23 (0.05) 15.2 (1.9) 140 (5) 7.5 (4.7) 306 (32) 91 (57) US-Bar 12 869 (113) 185 (14) 115 (12) 116 (12) 245 (17) 181 (67) 69 (13) 23 (0.06) 12.1 (1.3) 126 (16) 8.5 (8.7) 286 (58) 226 (68) US-UMB 11 973 (132) 176 (8) 103 (13) 124 (9) 256 (11) 112 (36) 42 (31) 24 (10) 13.5 (1.5) 140 (9) 25.3 (11.2) 253 (32) 133 (68) US-Ha1 15 927 (215) 203 (12) 101 (16) 109 (10) 245 (20) 85 (45) 43 (30) 17 (0.05) 13.0 (1.8) 124 (12) 10.7 (10.8) 273 (51) 194 (110) CA-Cbo 7 1449 (276) 193 (17) 103 (13) 113 (6) 253 (14) 146 (49) 73 (18) 21 (0.05) 13.5 (1.6) 133 (22) 16.4 (9.7) 243 (40) 161 (100) CA-TPD 6 949 (110) 191 (5) 99 (7) 121 (7) 239 (13) 44 (24) 35 (30) 11 (0.04) 15.2 (0.9) 132 (9) 22.2 (13.7) 255 (29) 231 (52) US-NR1 ENF 16 561 (69) 207 (12) 169 (8) 108 (10) 254 (16) 367 (63) 110 (18) 51 (13) 3.5 (1.2) 146 (11) 49.9 12.9 321 (29) 43 (30) US-Ho2 13 981 (161) 266 (16) 102 (13) 78 (10) 264 (16) 114 (57) 44 (32) 19 (0.07) 7.9 (1.1) 111 (14) 8.4 (8.2) 180 (37) 265 (101) US-Ho1 15 883 (93) 255 (14) 104 (14) 86 (10) 264 (9) 127 (56) 55 (28) 21 (0.07) 9.6 (0.9) 111 (13) 9.3 (7.6) 186 (17) 247 (97) US-Vcm 4 523 (100) 247 (16) 132 (15) 80 (12) 270 (6) 175 (77) 43 (25) 36 (0.08) 5.1 (2.3) NA (NA) 52.6 (11.1) 364 (24) 35 (9) US-Ha2 15 836 (146) 271 (17) 98 (17) 76 (12) 256 (20) 87 (45) 42 (30) 17 (0.06) 9.1 (1.3) 105 (22) 9.1 (8.3) 203 (42) 252 (85) CA-TP3 7 1301 (159) 269 (19) 99 (7) 83 (11) 285 (32) 49 (26) 38 (26) 12 (0.04) 8.8 (2.0) 92 (13) 16.8 (8.3) 157 (59) 203 (105) CA-TP4 7 979 (175) 267 (20) 99 (7) 85 (10) 256 (23) 48 (25) 38 (26) 12 (0.04) 11.1 (1.9) 93 (10) 16.4 (8.4) 213 (48) 313 (95)