Analyzing the Effects of Growing Season Length on the Net Ecosystem Production of an Alpine Grassland Using Model–Data Fusion

Alpine ecosystems are, similar to arctic ecosystems, characterized by a very long snow season. Previous studies investigating arctic or alpine ecosystems have shown that winter CO2 effluxes can dominate the annual balance and that the timing and duration of the snow cover plays a crucial role for plant growth and phenology and might also influence the growing season ecosystem CO2 strength and dynamics. The objective of this study was to analyze seasonal and annual CO2 balances of a grassland site at an elevation of 2440 m a.s.l in the Swiss central Alps. We continuously measured the NEP using the eddy covariance method from June 2013 to October 2014, covering two growing seasons and one winter. We analyzed the influence of snow melt date on the CO2 exchange dynamics at this site, because snow melt differed about 24 days between the 2 years. To this end, we employed a process-based ecosystem carbon cycling model to disentangle the co-occurring effects of growing season length, environmental conditions during the growing season, and physiological/structural properties of the canopy on the ecosystem carbon balance. During the measurement period, the site was a net sink for CO2 although winter efflux contributed significantly to the total balance. The cumulative growing season NEP as well as mean and maximum daily CO2 uptake rates was lower during the year with the later snow melt, and the results indicated that the differences were mainly due to differing growing season lengths.


INTRODUCTION
Long-term observations provide evidence that the earth's climate is changing. For example, an increase in global average temperature, more frequent and intense extreme climate events, and decreases in snowpack and snow cover have been documented (IPCC 2013), although locally these trends can be reversed especially at high latitudes or elevation (for example, Beniston and others 2003). According to the IPCC, these changes can be attributed to increased anthropogenic emissions of greenhouse gases (GHG), carbon dioxide (CO 2 ) being the most important of these GHGs.
Terrestrial ecosystems and the climate system are closely coupled through the carbon (C) cycle. The IPCC (2013) estimates that an increase in the C land sink has accumulated about 30% of the anthropogenic CO 2 emissions, reducing increases in the atmosphere by about 75 ppm between 1750 and 2011. To predict future changes in atmospheric CO 2 and global climate, it is crucial to quantify net ecosystem exchange of CO 2 (NEE) for different types of ecosystem and to understand how NEE is affected by CO 2 and climate.
Grassland ecosystems cover about 24% of the terrestrial area (Olson and others 1983) with different climatic conditions and management practices, and their role in the C cycle cannot be generalized. Measurements from various sites have shown that they can act as sources or sinks for CO 2 (Gilmanov and others 2010). Within individual sites, high interannual variability has also been observed and uncertainties remain as to the factors controlling the C dynamics in grassland ecosystems. Moreover, at least within Europe, there is a lack of data on year-round C fluxes of grassland sites at very high elevation.
High alpine ecosystems are-similar to ecosystems at high latitudes-characterized by a very long snow season. Although there is generally no or very low vascular plant growth during this period, microbial activity can continue under the snow pack (Gavazov and others 2017;Monson and others 2006;Panikov and others 2006) and for arctic environments it has been shown that winter CO 2 effluxes can dominate the annual CO 2 budget (Aurela and others 2002;Commane and others 2017;Euskirchen and others 2012).
The timing and duration of the snow cover affect plant growth and performance and, in the long term, also influence vegetation composition and sites with a long persistence of snow cover often support plant communities with lower growth rates and plant cover (Galen and Stanton 1993;Jonas and others 2008;Pickering and others 2014;Walker and others 1994). Spatial variation in snow cover can lead to productivity gradients or a mosaic of plant communities even within a small area when a typical snow cover pattern persists over time. Short-term responses to changes in the snow cover are complex and vary geographically and among plant communities (Humphreys and Lafleur 2011;Wipf and Rixen 2010). For many alpine and arctic plant species the timing of phenological events is linked to the timing of snow melt with earlier snow melt leading to, e.g., advanced flowering (Dunne and others 2003) and higher seed mass because of more time available for seed growth (Galen and Stanton 1993), whereas delayed snow melt can lead to delayed vegetation onset, delayed phenological phases, reduced biomass production and reproductive success of some vascular plants (Cooper and others 2011;Vitasse andothers 2017, 2010) but also to increased productivity in forbs (Wipf and Rixen 2010). Negative impacts of earlier snow melt on growth and productivity were also observed and linked to exposure to very cold temperatures without the insulating effect of snow cover or generally unfavorable spring environmental conditions (lower incident radiation or temperatures) when plant development starts (Bokhorst and others 2011;Galvagno and others 2013;Wipf and others 2009). A few studies also focused on the effect of snow melt timing on the C balance of alpine and arctic ecosystems. Although for some sites an earlier snow melt led to an increase in the annual net CO 2 uptake (Aurela and others 2004;Rennermalm and others 2005), no or even a reverse trend between net ecosystem production (NEP) and date of snow melt was found for others (Galvagno and others 2013;Humphreys and Lafleur 2011;Parmentier and others 2011). Although an earlier snow melt can reduce off-season emissions, Galvagno and others (2013) also observed that an early melt out date led to lower CO 2 uptake rates during the growing season, which was explained by structural and functional changes in the canopy (lower leaf area index (LAI) and chlorophyll content compared to years with a later snow melt). Humphreys and Lafleur (2011) also suggested that interannual differences in NEP are rather caused by interannual differences in canopy photosynthetic capacity (i.e., photosynthetic capacity derived at the canopy level, which is a function of LAI and leaf level photosynthetic capacity).
To analyze the effect of snow melt timing on the CO 2 balance of seasonally snow-covered ecosystems multiple interacting factors have to be disentangled. The snow cover duration dictates the timing/length for off-season emissions and the potential length of the growing season. It can influence the time available for plant development, but also determines the prevailing environmental conditions when plant development starts and, to a certain extent, during the entire growing season. Finally, the weather conditions during the growing season directly affect the ecosystem CO 2 balance, Effects of Growing Season Length on NEP which can interact with effects of the snow melt timing.
To master the multiple confounding factors influencing the CO 2 balance, process-based modeling approaches provide an important tool by providing complementary information to observational data, in particular on ecosystem attributes and processes that might not be directly measureable. Combining models and measured data through estimating model parameters based on the data by calibration, inverse modeling, or data assimilation techniques can improve model performance (Williams and others 2005) and the quantity and quality of the acquired information.
Here, we used the eddy covariance (EC) method to continuously measure the NEE at a grassland site in the Swiss central Alps at an elevation of 2440 m a.s.l. covering two growing seasons in 2013 and 2014. The objectives of this study were to establish annual and seasonal CO 2 balances and to analyze differences in the canopy photosynthetic capacity within the footprint of the measurements. In addition, we made use of a 'natural experiment' to investigate the impact of a delayed melt out on the net CO 2 uptake, because the site became snow free about 24 days later in the first study year compared to the second. To this end, we employed a processbased ecosystem C model based on the DALEC model (Fox and others 2009;Williams and others 2005), which was calibrated in a Bayesian framework.
We hypothesize (1) that on an annual basis, constrained by the short growing season, this high elevation grassland acts as a small sink for CO 2 at the most and (2) that a later spring snow melt decreases CO 2 uptake because of a shorter growing season and lower canopy photosynthetic capacity.

Site Description
The field site is an alpine pasture near Furka Pass (46°34¢ 36¢¢ N, 8°25¢ 17¢¢ E) in the Swiss central Alps at an elevation of about 2440 m a.s.l. on a southeast facing slope with an inclination of about 10-15°. Average annual precipitation and annual mean temperature at a close by location (Gü tsch, 17 km NE, 2283 m a.s.l.) were around 1450 mm and 0.4°C during the period 1981-2010. For the actual site, year-round meteorological data are only available since July 2012 and the mean precipitation during July-August from 2012 to 2016 was 270 mm, the mean air temperature during the same period was 8.3°C. Snow melt usually occurs in June and the growing season lasts about 2.5-3.5 months, with occasional snow events possible year round (Inauen and others 2012). The wind system is bimodal with westerly winds dominating during the night, and both, westerly and easterly winds, occurring during the day. The vegetation west of the EC tower consists of a Nardus community dominated by Nardus stricta and Carex curvula with some wetter patches of Eriophorum scheuchzeri and of Deschampsia cespitosa tussocks, whereas to the east the vegetation is sparser and typical snowbed vegetation dominated by Alchemilla pentaphyllea, Salix herbacea, Soldanella pusilla, Gnaphalium supinum, and Sibbaldia procumbens predominates. The soil was classified as partly podzolized alpine brown earth on siliceous bedrock (Inauen and others 2013).

CO 2 Flux Measurements
NEE was measured using the EC method others 2000, 2012;Baldocchi 2014). Measurements took place from June 21, 2013-October 8, 2014, covering two snow-free periods. The three orthogonal wind components and sonic temperature were measured using a 3-D sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA). CO 2 and H 2 O dry mole fractions were measured using an enclosed path infrared gas analyzer (IRGA) (LI-7200, LI-COR Inc., Lincoln, NE, USA) with an intake tube of 0.5 m length. The instruments were mounted at 3.5 m above ground, and data were sampled at 20 Hz on a data logger (CR5000, Campbell Scientific, Logan, UT, USA) and stored on a 2 GB data card. Because of lightning-caused instrument failure, the LI-7200 was temporarily replaced by an open-path gas analyzer (LI-7500, LI-COR Inc., Lincoln, NE, USA) during the period from August 19-October 8, 2013.
Fluxes were calculated according to commonly accepted procedures using the software EddyPro (LI-COR Inc., Lincoln, NE, USA). In brief, CO 2 fluxes were derived from the covariance of the vertical wind speed w and the CO 2 dry mole fraction c: F c ¼ w 0 c 0 , where primes denote fluctuations around the mean and the overbar a time average, in our case 30 min. Raw data procedures prior to flux calculation included de-spiking, a double rotation of wind data as well as the detection and compensation of time delays in the CO 2 signal. Lag times were estimated for each flux averaging period by maximizing the covariance between the CO 2 signal and the vertical wind component within a certain window of plausible time lags.
During winter time the site was not accessible for extended periods of time (because of high risk of avalanches in winter 2013/2014), making raw data storage on the data card impractical. Therefore, 30 min averages of the raw, unprocessed data of all variables as well as their standard deviations and covariances were computed and stored by the data logger from October 8, 2013-July 1, 2014, whereas during summer half-hourly data were logged as well as raw data. Covariances calculated from raw data after the different processing steps were compared to the unprocessed covariances. The postprocessing steps increased the covariance between the vertical wind speed and the CO 2 dry mole fraction by around 12%, mainly because of the compensation of the time delay of the CO 2 signal. Therefore, we multiplied the CO 2 fluxes calculated based on half-hourly covariances of unprocessed data (October 8, 2013-July 1, 2014) with a factor of 1.12.
Low-frequency losses due to the finite averaging time were corrected according to Moncrieff and others (2004). Correction for flux spectral losses in the high-frequency range due to the limited response of the instruments, sensor separation, path averaging, and so on, was done using the method by Fratini and others (2012). When no raw data were available, high-frequency losses were corrected according to Massman (2000).
Half-hourly flux data were discarded if the IRGA or sonic anemometer malfunctioned because of technical problems or weather conditions (precipitation, dew, snow). Moreover, time periods with the stationarity test or the deviation of the integral turbulence characteristics (Foken and Wichura 1996) exceeding 100% were excluded from analysis.
NEE was calculated as the sum of the vertical eddy flux and the storage term. The latter was estimated from the CO 2 concentrations based on a 1-point profile at the reference height (> 95% of all times less than ± 0.15 lmol m -2 s -1 ).
Nighttime NEE data were excluded when friction velocity (u * ) was below 0.18 m s -1 to avoid potential underestimation of ecosystem respiration during stable, calm nights when the assumptions underlying the eddy covariance method are not fulfilled. The u * threshold was estimated based on the Moving Point Test, a method for an objective threshold determination devised by Gu and others (2005).
The random uncertainty was estimated based on the neighboring days approach devised by Hollinger and Richardson (2005) by analyzing the difference of NEE values collected under similar conditions (same time of day, half-hourly photosynthetic photon flux density (PPFD) values within 100 lmol m -2 s -1 , air temperature (T a ) within 3°C, soil temperature (T s ) within 2°C, relative humidity (RH) within 20%, and wind speed within 1 m s -1 ) on two successive days.
For the analysis, net ecosystem exchange was expressed as net ecosystem production (NEP = -NEE), which equals the net gain or loss of C by the ecosystem (Reichle 1975), that is, positive values of NEP represent a net uptake of CO 2 by the ecosystem.
Overall, about 60% of the measured flux data were rejected because of quality criteria. This is a rather high rate, although it is quite common that 25-60% of eddy covariance data are omitted (for example, Papale and others 2006;Wohlfahrt and others 2008). To establish seasonal and annual budgets, data gaps were filled using common methodology. Data gaps less than 2 h were filled by linear interpolation. Longer gaps were filled by employing the functional relationship between NEP and soil temperature (T s , measured at 0.1 m depth) (Arrhenius-type relationship) for nighttime data (PPFD = 0) and between NEP and PPFD (Michaelis-Menten-type relationship, or so-called light response curves) for daytime data (PPFD > 0). For this purpose, the following models were fitted to biweekly blocks of half-hourly data shifted each 5 days: where a (lmol lmol -1 ) is the quantum yield, GPP sat (lmol m -2 s -1 ) the asymptotic value of the gross primary production (GPP) at high irradiance (photosynthetic capacity at canopy level), and R eco (lmol m -2 s -1 ) the ecosystem respiration, which is the sum of autotrophic (R a ) and heterotrophic (R h ) respiration. Because the vegetation differed in the east and west of the flux tower, a, GPP sat , and R eco and 95% confidence intervals were also estimated for the east (wind direction £ 180°) and west (wind direction > 180°) separately (for the analysis only and not used for gap filling). During the time with a permanent snow cover, no clear relationship between NEP and any environmental driver was found. Larger gaps were therefore filled with the mean value for the specific half-hour estimated from biweekly data blocks [mean diurnal variation; (Falge and others 2001)].
Based on the seasonal time-series of NEP we differentiated two phases: the carbon uptake period Effects of Growing Season Length on NEP (CUP) when the site represented a sink of CO 2 , defined as the period when the 5-day mean of NEP is positive (Galvagno and others 2013;Wohlfahrt and others 2013), and the period during which the site represented a source of CO 2 .

Ancillary Measurements
Meteorological conditions, including air temperature (T a ) and relative humidity (RH), global radiation (i.e., diffuse and direct solar radiation; R g ), precipitation (P) and snow height (SH), were monitored by a weather station set up at 2.4 m above ground. In addition, soil temperature (T s ) at 0.1 and 0.25 m below ground and soil moisture (soil water content (SWC) in vol%) at 0.05 and 0.15 m below ground were measured. If not indicated otherwise, T s at 0.1 m and SWC at 0.05 m were used. Because photosynthetic active radiation was not directly measured, PPFD (lmol m -2 s -1 ) was approximated as PPFD = 2 * R g (W m -2 ) (Gonzá lez and Calbó 2002). The fraction of diffuse radiation was estimated using an empirical relationship between R g and the clearness index parameterized based on data from a study site in the Eastern Alps (Wohlfahrt and others 2016).
Because the area under the snow height sensor was among the first to become snow free, while the rest of the field site was still covered by snow, snow cover duration was determined from webcam images.

Models
We used a box model based on the Data Assimilation Linked Ecosystem Carbon (DALEC) model (Fox and others 2009;Williams and others 2005) to represent the C cycle of the alpine grassland. The model consists of six C pools that represent the C content of foliage (C f ), roots (C r ), necromass (standing dead biomass) (C n ), litter (C lit ), and soil organic matter (C som ). At our grassland site, many of the dominant species are hemicryptophytes and therefore the model contains a labile C pool (C lab ) that represents the overwintering buds and supports leaf formation in spring. The pools are connected via fluxes as illustrated in Figure 1, where A denotes allocation to a pool, L the litterfall from the respective pool, and D the decomposition of litter. It was assumed [similar to Williams and others (2005)], that autotrophic respiration is a fixed fraction of GPP and that all C fixed during the day is either respired (autotrophic respiration) or allocated to one of the plant tissue pools (C f , C r , or C lab ), there is no direct environmental influence on allocation/litterfall and those fluxes are donorcontrolled with constant rate parameters, soil transformations are temperature-sensitive (see Appendix A) and there is no C loss due to herbivory or leaching.
Preliminary data analysis and model runs indicated that leaf senescence was predominantly driven by photoperiod, as was found in previous studies (Kö rner 2003). Therefore, in the model a minimum day length threshold, which was estimated during model calibration, controlled leaf senescence in fall.
For a detailed description of the model equations see Appendix A.
The daily C input into the system, i.e., GPP, was calculated as a function of the green plant area index (GAI, area of green plant matter per ground area; directly estimated from C f ) by use of a sun/ shade big-leaf model of canopy photosynthesis based largely on de Pury and Farquhar (1997), with modifications according to Goudriaan and Laar (1994), Wang and Leuning (1998) and Smith (1937) The model contained 17 unknown parameters, and the starting values of the six C pools were also unknown. Three of the parameters associated with the GPP model were fixed to values taken from the literature (Table 1), the remaining 14 parameters and initial conditions were estimated using DREAM (Vrugt and others 2009) to reconcile model output with measured NEP. For this purpose, two data sets of daily values-NEP day and NEP night -were generated from the non-gap-filled data; daily values were calculated for days when at least 80 and 50% of the half-hourly NEP measurements during day and night were valid, respectively.
DREAM is a multi-chain Markov Chain Monte Carlo (MCMC) algorithm for statistical inference of parameters using Bayesian statistics. A Bayesian calibration estimates the joint probability distribution of all parameters conditional on the available data rather than estimating specific values for each parameter. This so-called posterior distribution is computed from the prior distribution of the parameters and their respective likelihood given the observed data. The prior distribution is determined by the information on the parameters prior to data collection and analysis. The likelihood function provides a measure of how well the model fits the data. For most complex process-based models the posterior distribution cannot be determined analytically and to estimate the likelihood the model needs to be run. This can be realized using a sampling method like MCMC, which generates a random walk through the parameter search space. The basis of this method consists of the following steps: first, based on the starting point, a candidate parameter set is sampled from a proposed distribution. Then the likelihood is estimated and the candidate is either accepted or rejected according to an acceptance probability. If the candidate is accepted, the chain moves to the new location; otherwise, it remains in the current location. Repeating those steps eventually results in a Markov chain representing the target distribution. To explore multi-dimensional parameter spaces rapidly and adequately, multiple chains are run in parallel and achievement of convergence to a limiting distribution is estimated (Gelman and Rubin 1992).
In our case, a uniform prior distribution of the parameters was used with lower/upper bounds (given in Table 1) determined based on information from the literature and initial test runs of the model.
To account for non-normality, heteroscedasticity and correlation of model residuals, a common problem in ecological modeling, a generalized likelihood function (Schoups and Vrugt 2010) was applied and the appropriate residual error distribution was determined concurrently with the model parameters. To estimate model uncertainty due to parameter uncertainty, model output was generated with 1000 parameter sets from the posterior distribution and the 97.5 and 2.5% prediction percentiles were calculated. Total predictive uncertainty was estimated by adding the modeled residuals based on the estimated residual error distribution. Collinearity of the model parameters was checked by computing the correlation matrix of the parameter estimates after the chains converged to the target distribution.
To investigate the differences between the two study years, the model was calibrated against NEP data from the individual years, resulting in two sets of model parameters ('para13' and 'para14'). However, the initial condition of the C lab , C r , C lit , and C som pools was set equal for both years, where the initial magnitude was estimated by calibrating the model against the entire NEP data set of both years. The initial magnitudes of C f and C n were set to zero, because simulations commenced before onset of the growing season.
Subsequently, four forward model runs were performed, using the two parameter sets and the environmental conditions of the 2 years in a factorial design (Table 2).
Finally, to analyze the direct effects of growing season weather during 2013 and 2014 on NEP without the complication of seasonally changing C pool sizes, the model was run using each year's weather conditions and the 2013 parameters, but instead of simulating the dynamics of the C pools, their values were set equal to the C-pool values of the 2013_para13-model output. Thus, the development of the C pools and therefore also GAI was forced to be identical during both years and differences in NEP between the 2 years were only due to differences in environmental conditions (incident radiation and temperature). Simulated values for NEP, GPP, and R eco were analyzed during a period that was snow free in both years and also covered the peak uptake phase of the 2 years (July 15-August 15).   Assumptions for parameters that were fixed during the calibration are given in column 4; for parameters that were optimized, lower/upper bounds of prior distribution are given in column 4. a NPP (net primary production) = GPP-R a .

Environmental Conditions
The seasonal course of incoming global radiation showed some lower values especially during the summer months in 2014 compared to 2013 (Figure 2A), which was due to more foggy and overcast weather conditions. The mean temperature during July-August was 8.6°C in 2013 and 6.1°C in 2014. Highest temperatures were reached around midend of July, well after the summer solstice. July-August precipitation was 300 and 390 mm in 2013 and 2014, respectively. SWC was slightly lower during the main growing season in 2013. However, SWC was never below 17% during the growing season in both years ( Figure 2D), whereas field capacity was estimated to around 22% indicating that there was no water limitation (Boyer and Kramer 1995). T s stayed around zero during the snowcovered period, but rapidly increased after snow melt and followed the general course of T a afterward. In both years, snow height reached maximum values of about 1.30 m ( Figure 2E). In 2013, the snow height started to continuously decrease in the beginning of June and on June 19 the area directly under the snow height sensor was snow free. In 2014, snow melt already started in the middle of May and the first snow-free day was registered by the sensor on June 8. The day when a major part of the grassland was snow free to a similar extent in 2013 and 2014 determined using webcam images was July 3, 2013, and June 9, 2014. Thus, there was a difference of 24 days in snow melt date between the 2 years.

Net Ecosystem Productivity
Because the vegetation appeared sparser in the east compared to the west of the measuring tower, NEP in relation to wind direction was analyzed. During the measurement period, the wind was mainly blowing from the east or west, with CO 2 uptake rates being slightly higher in the west ( Figure C1 in Appendix C). Also, the estimated values for GPP sat were higher in the west, especially in the beginning of the growing season 2013, whereas the estimated values for R eco hardly showed a difference (Figure C2). Comparing the 2 years, GPP sat and R eco were higher in 2014.
Because the difference between east and west appeared small, measurements from all wind directions were pooled for the rest of the analysis.
During the entire measurement period (i.e., two growing seasons and one winter), the cumulative NEP was about 160 g C m -2 (Figure 3). Considering one full year with data coverage, i.e., from July 2013 to July 2014, the site represented a small net CO 2 sink (20 g C m -2 ).
In both years, the CUP started 14 days after snow melt and lasted from July 17, 2013-September 26, 2013 (72 days), and from June 23, 2014-September 28, 2014 (98 days). During those periods, about 100 and 150 g C m -2 were taken up with an average uptake rate of 1.40 and 1.56 g C m -2 d -1 , respectively. The maximum daily NEP of 3.16 g C m -2 d -1 in 2013 was also lower than in 2014 (4.11 g C m -2 d -1 ). The later start of the CUP in 2013 is reflected in the monthly balance, where the cumulative NEP in July was considerably lower in 2013 than in 2014, whereas the difference between the 2 years was much lower in the following 2 months (Figure 3).
At the time of snow melt the NEP was about the same during both years ( Figure 4A). Two weeks after the snow melt date, the mean diurnal variation of NEP showed a clear pattern of CO 2 uptake during the day and CO 2 release during the night, with slightly higher uptake in 2014 ( Figure 4B). During the peak uptake, which was reached 34 days after snow melt in 2013 and 47 days after snow melt in 2014, both CO 2 uptake and respiratory loss were higher in 2014 ( Figure 4C), despite the mean diurnal course of PPFD being lower in 2014 ( Figure 4F).
The light response curves also showed higher CO 2 uptake at the beginning of the CUP as well as during the peak uptake phase in 2014 ( Figure C3). Table 2. Four possible forward model runs using the environmental conditions and the parameter sets of the two individual years in a factorial design During the non-CUP between the 2 years, about 90 g C m -2 was released with a mean daily rate of -0.32 g C m -2 d -1 . Highest CO 2 release was observed after the end of the CUP before the permanent snow cover established and directly after snow melt with a maximum of -1.49 g C m -2 d -1 after the end of the CUP 2013. After establishment of the snow cover, daily CO 2 release rates were relatively stable with an overall declining trend over time ( Figure 5).

Simulation Analyses
To investigate whether the lower C uptake in 2013 was mainly due to the shorter growing season, an effect of biotic response to spring environmental conditions, or the direct effect of the weather conditions during the growing season, the process- Figure 2. Annual course of A daily sum of global radiation; B daily mean air temperature; C daily mean soil temperature; D daily soil water content; and E daily snow height for the two study years. based ecosystem C model was utilized. Model performance was very good during daytime and satisfactory during nighttime (Nash-Sutcliffe efficiency (NSE) ‡ 0.9 and ‡ 0.5, respectively (Moriasi and others 2007); Figure 6). For 2013, the mean absolute error (MAE) between daytime model and data values was 0.29 g C m -2 d -1 (slope of 0.85, intercept of 0.03 g C m -2 d -1 ; R 2 = 0.92; Figure 6B). For 2014, the MAE was 0.30 g C m -2 d -1 (slope of 0.84, intercept of 0.08 g C m -2 d -1 ; R 2 = 0.92; Figure 6F). For nighttime values, the MAE was 0.19 g C m -2 d -1 (slope of 0.6, intercept of -0.16 g C m -2 d -1 ; R 2 = 0.52) for 2013 (Figure 6D) and 0.25 g C m -2 d -1 (slope of 0.6, intercept of -0.15 g C m -2 d -1 ; R 2 = 0.65) for 2014 ( Figure 6H). The estimates for the initial conditions of the C pools as well as the optimized model parameters are listed in Table 1.
The correlation analysis indicated high collinearity between some of the model parameters ( Figure 7). For example, A max was negatively correlated with the allocation rate of NPP to C f (r = -0.8), as was the turnover rate of roots with the mineralization (r = -0.9) and decomposition (r = -0.7) rates of litter. The allocation rate to C n was negatively correlated with the allocation rate to C lab (r = -0.7) and the turnover rate of C n (r = -0.7). The fraction of GPP respired was positively correlated with the fraction of C lab allocated to C f (r = 0.7), as was the allocation rate to C lab with the turnover rate of C n (r = 0.6). The mineralization rate of litter was positively correlated with the decomposition rate of litter (r = 0.7).
Most parameters were in the same order of magnitude for the 2 years. However, the allocation rate to C f as well as the turnover rate of C f in 2013 was about twice as large as in 2014. On the other hand, the turnover rate of C r was considerably smaller in 2013 compared to 2014. In fall, the model parameters showed a lower allocation to C lab in 2013, while the allocation to C n was higher compared to 2014. The allocation rate from C lab to C f , and therefore spring start-up, was about three times larger in 2014 than in 2013. This resulted in a slightly steeper increase in the simulated GAI during the beginning of the growing season in 2014 compared to 2013. However, the maximum GAI was about the same during both years (2.3 m 2 m -2 ) and was reached on August 19, 2013, and on August 9, 2014. These results also compare well with GAI data from MODIS satellite observations (ORNL DAAC 2008) of the area around the flux tower (1 km resolution) (Figure 8).
For all model runs, GAI increased more rapidly at the beginning of the growing season with 2014 parameters ( Figure C4). However, under the same environmental conditions, the 2013 parameters resulted in higher maximum GAI, with the highest value for 2014 environmental conditions and 2013 parameters (2014_para13-simulation).
This pattern was also reflected in the simulated NEP (Figure 9). In both years, using the 2014 parameters led to a slightly steeper increase in NEP at the beginning of the CUP, while at the end of the CUP the cumulative NEP was higher with 2013 parameters when comparing the model output of the runs under the same environmental conditions. The maximum daily NEP was 3.2, 2.8, 5.0, and 3.5 g C m -2 d -1 for the 2013_para13-, 2013_ para14-, 2014_para13-, and 2014_para14-model runs, respectively. The annual cumulative NEP was about the same for both parameter sets within the Because NEP is simulated as a function of GAI, C lit , C som , PPFD, and T s , the results above reflect the combined effect of ecosystem status and actual weather conditions (i.e., differences in peak uptake could result from both, a difference in ecosystem status (e.g., GAI, photosynthetic capacity) or a difference in for example incoming radiation). The analysis of the direct impact of growing season weather on NEP, however, showed no effect (data not shown).

DISCUSSION
In high mountain ecosystems, plant development is very dependent on abiotic factors, for example, the snow cover plays a crucial role in determining the potential growing season (Galvagno and others 2013;Inouye and Wielgolaski 2013;Kö rner 2003;Wohlfahrt and others 2013). In this study, we used EC measurements and a modeling approach to investigate the magnitude and differences in the NEP of an alpine grassland during two growing seasons differing 24 days in snow melt date. During the snow-free season, the grassland was a net sink for CO 2 during both years with total uptake com-parable to sites at very high elevation, for example, the Tibetan Plateau (Kato and others 2006), or high latitude (Humphreys and Lafleur 2011). The time between snow melt and the start of the CUP was the same during both years (14 days) indicating that the onset of vegetation is tightly linked to snow melt at this site. Despite differences in snow melt date, peak and end of the CUP were reached around the same day of the year during both years, thus leading to a shorter CUP in the year with delayed snow melt (2013).
Cumulative NEP during the CUP as well as mean and maximum daily NEP rates were all higher for the year with the longer CUP. In contrast to our results, Galvagno and others (2013) found that an earlier snow melt led to no change in the cumulative growing season CO 2 uptake despite a longer CUP because of lower CO 2 uptake rates during the elongated CUP. However, they investigated a grassland at lower elevation where the snow melt date occurred well before the summer solstice. In that case, an earlier snow melt date means a shift toward shorter day length, which has also been found to cause a lengthening of the period between snow melt and the start of the CUP (Wohlfahrt and others 2013).
At high elevation as in our case, the snow melt usually occurs around the summer solstice and the determination of the growing season length by the timing of the snow melt in spring seemed to be the dominant factor for the growing season ecosystem CO 2 balance.
The end of the CUP was reached well before a permanent snow cover established and senescence appeared to be driven predominantly by photoperiod (peak and end of CUP were reached around the same day of the year in both years despite differences in temperatures). Highest CO 2 emission rates were observed directly after the end of the CUP and also relative high rates directly after snow  melt before the start of the new CUP. This result is in line with findings by Aurela and others (2002) who observed maximum daily net effluxes before and after the sink period in a subarctic fen. Commane and others (2017) noted the importance of winter respiration especially during early winter when soil temperature is still relatively high in the arctic.
During the snow season, CO 2 release continued with declining CO 2 emission rates over time. Because of the snow cover, the soil temperature was decoupled from air temperature and continuously stayed around -0.2 to -0.1°C. Thus, the data indicate that microbial activity continued under the snow cover and CO 2 fluxes were determined by substrate limitation (Liptzin and others 2009). Although the average winter flux of -0.33 g C m -2 d -1 was rather low, winter efflux contributed significantly to the annual budget because of the long snow season at this high elevation.
The modeling results indicated that the site was a net sink on an annual basis, with a cumulative NEP of 35 g C m -2 in 2013 and 130 g C m -2 in 2014. Although the model performed well in simulating daytime NEP, comparing simulated and measured nighttime fluxes showed relatively large MAE compared to the observed data range. This result could indicate that the chosen model structure does not fully capture relevant ecosystem processes in a way that parameters can be optimized for daytime and nighttime fluxes simultaneously, which could be improved in future studies. However, when using the eddy covariance method, nighttime flux measurements can also be expected to have higher uncertainty (Aubinet and others 2012) and contributing to higher discrepancies between simulated and measured values. Both, simulated nighttime and daytime data tended to be underes-timated. Correcting for this bias resulted in the site being a slightly larger net sink annually (65 g C m -2 in 2013 and 135 g C m -2 in 2014).
The larger sink in 2014 compared to 2013 resulted from both, a shorter snow season with net CO 2 efflux as well as the longer CUP with a slightly higher average CO 2 uptake rate.
Comparing the simulated NEP of the four model runs showed that within the same year, that is the same meteorological forcing, the annual NEP was about the same for both parameter sets. However, comparing the 2 years, the NEP was always lower under 2013 environmental conditions, which were characterized by the later snow melt and drier and warmer conditions during the growing season. GPP and R eco , and therefore NEP, are processes that depend on weather conditions. Our analysis of the direct effect of growing season weather on NEP did not show any significant differences between the two study years. Thus, the modeling results reinforce that the lower cumulative CO 2 uptake during the year with a longer persistence of the snow cover in spring was primarily due to the shortened growing season.
The simulated GAI also indicated differences in the development of the canopy between the 2 years. The earlier snow melt in 2014 enabled a more rapid development of the grassland (i.e., faster increase in GAI) at the beginning of the growing season. However, a faster rate during the main growing season in 2013 led to the maximum GAI being about the same in both years. Applying the 2013 model parameters under 2014 environmental conditions (2014_para13-simulation) indicated that with an elongated growing season the GAI would have reached a much higher value. However, a GAI of almost 4 is unrealistic at this site and is due to not considering other potentially growth GAI directly influences GPP, because it is a measure of the photosynthetically active plant tissue. The reason why the simulated annual NEP of the 2014_para13-and 2014_para14-runs was about the same, despite the maximum GAI being considerably higher in the 2014_para13-simulation, is twofold. On the one hand, respiratory losses were higher for the para13-run. On the other hand, GPP depends not only on the GAI, but also on the leaf photosynthetic capacity, which was estimated to be lower in 2013. This pattern, however, could also be a result of parameter equifinality-a common problem in ecological models due to overparameterization and poor identifiability of model parameters (Medlyn and others 2005). The results of our collinearity test showed high correlation between several parameters. For example A max and the rate of NPP allocated to the foliage pool, which directly influences the GAI, were highly correlated. Because no direct measurements of any of the C pools, GAI, or photosynthetic capacity were available, model parameters were estimated by fitting the model only to NEP data and for example multiple combinations of values for A max and the allocation rate to the foliage pool could potentially lead to an equally good fit. This highlights the fact that the results for those parameters cannot be interpreted unambiguously and stresses the importance of having multiple data sets to improve model output and reduce parameter uncertainty, as was also shown by Williams and others (2005). The results of the parameter optimization and parameter collinearity test can be used to determine which additional measurements are most valuable to improve modeling results and interpretation of those results for future studies.
Nevertheless, comparison of simulated C pools with values from the literature showed good agreement. C content of alpine grassland soils was analyzed by Budge and others (2011) at close by locations, and values of 1800-3150 g C m -2 in the upper soil layer (0-5 cm) were found with an estimated proportion of labile C between 71 and 86%. The C som starting pool of our model, which primarily represents the labile C content, was estimated as 1500 g C m -2 . Dry mass of litter was also measured by Budge and others (2011) and was between 800 and 2900 g m -2 . Assuming 50% C per dry matter, this measurement is higher than the estimated C litter starting pool of 100 g C m -2 . However, the measurements by Budge and others took place in fall when litter can be expected to be at maximum and their values for litter also included fine roots.
Above-ground peak biomass at our study site is estimated to around 200 g dry matter m -2 (E. Hiltbrunner, personal communication). Data on below ground biomass are sparse and very uncertain but often assumed to be in the range of or higher than above-ground biomass, which was for example also observed for some glacier forefield species near our study site (Inauen and others 2012). This observation falls well in line with the estimated C r (below ground biomass) starting pool of 120 g C m -2 and maximum simulated C f (aboveground biomass) pool of around 100 g C m -2 . Inauen and others (2013) measured peak LAIs of 0.8-2.6 m 2 m -2 for the most abundant vegetation types in the same study area in 2010, and peak values are generally reached around the beginning of August. Our simulated GAI reached peak values of 2.3 m 2 m -2 around mid-August.

CONCLUSION
Despite the very short growing season, the investigated alpine grassland was a net sink for CO 2 on an annual basis because of fast onset and increase in photosynthetic activity after snow melt. The combined analysis of NEP measurements and ecosystem-based model simulations showed that the snow melt date had a major effect on the seasonal and annual CO 2 balance, where the determination of the length of the growing season by the snow melt date in spring appeared to be the dominant factor. Thus, ecosystems at this elevation may benefit from an earlier snow melt leading to higher C uptake during a longer growing season. This gain, however, can be overshadowed when a major advancement of the snow melt date leads to a shift of the potential start of the growing season toward a period with shorter days or exposure of the vegetation to low temperatures.

OP E N A C C ES S
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/ 4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.