Growth Response of Grasslands to Snow Cover Duration Printer-friendly Version Interactive Discussion Growth Response of Temperate Mountain Grasslands to Inter-annual Variations of Snow Cover Duration Growth Response of Grasslands to Snow Cover Duration Printer-friendly Version Interactive Discussion

Introduction Conclusions References Tables Figures Back Close Full Screen / Esc This discussion paper is/has been under review for the journal Biogeosciences (BG). Please refer to the corresponding final paper in BG if available. Abstract Introduction Conclusions References Tables Figures Back Close Full Screen / Esc Abstract A remote sensing approach is used to examine the direct and indirect effects of snow cover duration and weather conditions on the growth response of mountain grasslands located above the tree line in the French Alps. Time-integrated normalized difference vegetation index (NDVI int), used as a surrogate for aboveground primary 5 productivity, and snow cover duration were derived from a 13 year long time series of the Moderate Resolution Imaging Spectro-radiometer (MODIS). A regional-scale meteorological forcing that accounted for topographical effects was provided by the SAFRAN–Crocus–MEPRA model chain. A hierarchical path analysis was developed to analyze the multivariate causal relationships between forcing variables and proxies 10 of primary productivity. Inter-annual variations in primary productivity were primarily governed by year-to-year variations in the length of the snow-free period and to a much lesser extent by temperature and precipitation during the growing season. A prolonged snow cover reduces the number and magnitude of frost events during the initial growth period but this has a negligeable impact on NDVI int as compared to the strong negative 15 effect of a delayed snow melting. The maximum NDVI slightly responded to increased summer precipitation and temperature but the impact on productivity was weak. The period spanning from peak standing biomass to the first snowfall accounted for two thirds of NDVI int and this explained the high sensitivity of NDVI int to autumn temperature and autumn rainfall that control the timing of the first snowfall. The ability of mountain 20 plants to maintain green tissues during the whole snow-free period along with the relatively low responsiveness of peak standing biomass to summer meteorological conditions led to the conclusion that the length of the snow-free period is the primary driver of the inter-annual variations in primary productivity of mountain grasslands. Abstract Introduction Conclusions References Tables Figures Back Close Full Screen / Esc


Introduction
Temperate mountain grasslands are seasonally snow-covered ecosystems that have to cope with a limited period of growth (Körner, 1999).The extent to which the length of the snowfree period controls the primary production of mountain grasslands is still debated.On the one hand, snow cover manipulation experiments and time series analyses of groundbased measurements have generally shown a decrease in biomass production under shortened growing season length (Wipf and Rixen, 2010;Rammig et al., 2010).On the other hand, several studies have pointed to the increasing risk of spring frost damage and summer water shortage following an early snowmelt and the associated detrimental effects on biomass production (Baptist et al., 2010;Ernakovich et al., 2014;Inouye, 2000).In addition, both soil microbial nitrogen immobilization and accumulation of inorganic nitrogen are enhanced under deep and long-lasting snowpacks (Brooks et al., 1998), and plants may benefit from increased flush of nutrients and ameliorated soil water balance following unusually long winters.To better understand the growth response of alpine grasslands to changing snow cover duration, it thus seems pivotal (i) to assess the contribution of the different components of the growth response, particularly the duration of the favorable period of growth and the peak standing biomass; (ii) to account for the effect of meteorologi-P.Choler: Growth response of grasslands to snow cover duration cal forcing variables on both snow cover dynamics and on plant growth, and (iii) to disentangle the direct and indirect effects, i.e., effects mediated by other forcing variables, of snow cover on land surface phenology and primary productivity.
From a phenomenological point of view, annual primary production may be viewed as the outcome of two things, namely the time available for biomass production and the amount of biomass produced per unit of time.For seasonally snow-covered ecosystems, this translates into two fundamental questions: (i) to what extent does the length of the snow-free period determine the length of plant activity and (ii) what are the main drivers controlling the instantaneous primary production rate of grasslands during the snow-free period?A number of studies have provided evidence for the non-independence of these two facets of growth response by noting that the biomass production rate increases when snow melting is delayed and that grasslands are able to partially recover the time lost when the winter was atypically long (Walker et al., 1994;Jonas et al., 2008).However, most of these studies focused on the initial period of growth -i.e., from the onset of greenness to the time of peak standing biomass -and therefore little is known about the overall relationship between the mean production rate and the total length of the snow-free period.Eddy covariance measurements have shown that the amount of carbon fixed from the peak standing biomass to the first snowfall represents a significant contribution to the gross primary productivity (GPP) (e.g., Rossini et al., 2012).Accounting for the full period of plant activity when examining how primary production of grasslands adjusts to inter-annual variations in meteorological conditions seems to be thus essential.
Remote sensing provides invaluable data for tracking ecosystem phenology over a broad spatial scale as well as inter-annual variations in phenological stages over extended time periods (Pettorelli et al., 2005).For temperature-limited ecosystems, numerous studies focused on arctic areas have established that the observed decadal trend toward an earlier snowmelt has translated into an extended growing season and enhanced greenness (Myneni et al., 1997;Jia et al., 2003).By contrast, the phenology of high elevation grasslands has not received the same degree of attention, partly because there are a number of methodological problems in using remote sensing data in topographically complex terrain, including scale mismatches, geolocation errors, and vegetation heterogeneity (Fontana et al., 2009;Tan et al., 2006).That said, some studies have used moderate-resolution imagery to document the contrasting responses of low and high vegetation to the 2003 heatwave in the Alps (Jolly, 2005;Reichstein et al., 2007) or to characterize the land surface phenology of high elevation areas in the Rockies (Dunn and de Beurs, 2011), the Alps (Fontana et al., 2008) or the Tibetan Plateau (Li et al., 2007).However, none of these studies have comprehensively examined the direct and indirect effect of meteorological forcing variables and snow cover duration on the different components of annual biomass production in mountain grasslands.
In this paper, I used remotely sensed time series of the Normalized Difference Snow index (NDSI) and of the Normalized Difference Vegetation Index (NDVI) to characterize snow cover dynamics and growth response of mountain grasslands.Time-integrated NDVI (NDVI int ) and the product of NDVI and photosynthetically active radiation (PAR) were taken as surrogates of aboveground primary productivity, while maximum NDVI (NDVI max ) was used as an indicator of growth responsiveness to weather conditions during the summer.My main aim is to decipher the interplay of snow cover dynamics, weather conditions and growth responsiveness affecting NDVI int .Specifically, I addressed three questions: (i) what is the relative contribution of the growing season length and NDVI max in determining the inter-annual variations in primary productivity?(ii) What are the direct and indirect effects of the snow cover dynamics on productivity?and (iii) What is the sensitivity of NDVI int to interannual variations in temperature and precipitation during the growing season?The study was based on 121 grasslandcovered high elevation sites located in the French Alps.Sites were chosen to enable a remote sensing characterization of their land surface phenology using the Moderateresolution Imaging Spectroradiometer (MODIS).Meteorological forcing was provided by the SAFRAN-CROCUS-MEPRA model chain that accounts for topographical effects (Durand et al., 2009c).I implemented a hierarchical path analysis to analyze the multivariate causal relationships between meteorological forcing, snow cover, and NDVIderived proxies of grassland phenology and primary productivity.

Selection of study sites
The selection of sites across the French Alps was made by combining several georeferenced databases and expert knowledge.My primary source of information was the 100 m-resolution CORINE land cover 2000 database produced by the European Topic Centre on Spatial Information and Analysis (European Environment Agency, 2007) that identifies 44 land cover classes based on the visual interpretation of high-resolution satellite images and from which I selected the 3.2.1 class corresponding to "Natural grasslands".Natural grasslands located between 2000 m and 2600 m above sea level were extracted using a 50 mresolution digital elevation model from the Institut Géographique National (IGN).I then calculated the perimeter (P ), area (A) and the mean slope of each resulting group of adjacent pixels, hereafter referred to as polygons, and kept only those that had an area greater than 20 ha, an index of compactness (C = 4π A/P 2 ) greater than 0.1, and a mean slope smaller than 10 • .The first two criteria ensured that polygons were large enough and sufficiently round-shaped to include several 250 m MODIS contiguous cells and to limit edge effects.The third criterion reduced the uncertainty in reflectance estimates associated with steep slopes and different aspects within the same polygon.Moreover, steep slopes usually exhibit sparser plant cover with low seasonal amplitude of NDVI, which reduces the signal-to-noise ratio of remote sensing data.Finally, I visually double-checked the land cover of all polygons by using 50 cm-resolution aerial photographs from 2008 or 2009.This last step was required to discard polygons located within ski resorts and possibly including patches of sown grasslands, and polygons too close to mountain lakes and including swampy vegetation.I also verified that all polygons were located above the treeline.

Climate data
Time series of temperature, precipitation and incoming shortwave radiation were estimated by the SAFRAN-CROCUS-MEPRA meteorological model developed by Météo-France for the French Alps.Details on input data, methodology, and validation of this model are provided in Durand et al. (2009a, b).To summarize, the model combines observed data from a network of weather stations and estimates from numerical weather forecasting models to provide hourly data of atmospheric parameters including air temperature, precipitation, and incoming solar radiation.Simulations are performed for 23 different massifs of the French Alps (Fig. 1), each of which is subdivided according to the following topographic classes: 300 m elevation bands, seven slope aspect classes (north, flat, east, southeast, south, southwest, and west) and two slope classes (20 or 40 • ).The delineation of massifs was based on both climatological homogeneity, especially precipitation, and physiographic features.To date, SAFRAN is the only operational product that accounts for topographic features in modeling meteorological land surface parameters for the different massifs of the French Alps.

MODIS data
The MOD09A1 and MOD09Q1 surface reflectance products corresponding to tile h18.v4 (40-50 • N, 0-15.6 • E) were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) (ftp://e4ftl01.cr.usgs.gov).A total of 499 scenes covering the period from 18 February 2000 to 27 December 2012 was acquired for further processing.Data are composite reflectance, i.e., represent the highest observed value over an 8-day period.Surface reflectance in the red (RED), green (GREEN), near-infrared (NIR), and midinfrared (MIR) were used to calculate an NDVI at 250 m following: NDVI = (NIR − RED)/(NIR + RED) (1) and an NDSI at 500 m using the algorithm implemented in Salomonson and Appel (2004): NDVI and NDSI values were averaged for each polygon.Missing or low-quality data were identified by examining quality assurance information contained in MOD09Q1 products and interpolated using cubic smoothing spline.NDVI or NDSI values that were 2 times larger or smaller than the average of the two preceding values, and the two following values were considered as outliers and discarded.Time series were gap-filled using cubic spline interpolation and smoothed using the Savitzky-Golay filter with a moving window of length n = 2 and a quadratic polynomial fitted to 2n + 1 points (Savitzky and Golay, 1964).
A high NDSI and low NDVI were indicative of wintertime whereas a low NDSI and a high NDVI were indicative of the growing season (Fig. 2).Here I used the criteria NDSI / NDVI < 1 to estimate the length of the snow-free period, hereafter referred to as P sf , at the polygon level (Fig. 2).This ratio was chosen as a simple and consistent way to set variables used in this study: date of snowmelt (TSNOW melt ), maximum NDVI (NDVI max ) and date of NDVI max (TNDVI max ), date of snowfall (TSNOW fall ), length of the snow-free period (P sf ), length of the initial growth period (P g ), length of the senescence period (P s ), and time-integrated NDVI over the growth period (NDVI int,g ) and over the senescence period (NDVI int,s ).
the start (TSNOW melt ) and the end (TSNOW fall ) of the snowfree period across polygons and years.Ground-based observations corresponding to one MOD09A1 pixel (Lautaret pass, 6.4170 • longitude and 45.0402 • latitude) and including visual inspection, analysis of images acquired with timelapse cameras and continuous monitoring of soil temperature and snow height, showed that this ratio provides a fair estimate of snow cover dynamics (Supplement Fig. S1).Further analyses also indicated that P sf is relatively insensitive to changes in the NDVI/NDSI thresholds with 95 % of the polygon × year combinations exhibiting less than 2 days of shortening when the threshold was set to 1.1 and less than 3 days of lengthening when the threshold was set to 0.9 (Fig. S2).Finally, changing the threshold within this range had no impact on the main results of the path analysis.The yearly maximum NDVI value (NDVI max ) was calculated as the average of the three highest daily consecutive values of NDVI and the corresponding middle date was noted TNDVI max .
The GPP of grasslands could be derived from remote sensing data following a framework originally published by Monteith (Monteith, 1977).In this approach, GPP is modeled as the product of the incident PA, the fraction of PAR absorbed by vegetation (fPAR), and a light-use efficiency parameter (LUE) that expresses the efficiency of light conversion to carbon fixation.It has been shown that fPAR can be linearly related to vegetation indices under a large combination of vegetation, soil-and atmospheric conditions (Myneni and Williams, 1994).Assuming that LUE was constant for a given polygon, I therefore approximated inter-annual variations in GPP using the time-integrated value of the product NDVI × PAR, hereafter referred to as GPP int , over the grow-ing season and calculated as follows: where T is the number of days for which NDVI was above NDVI thr .I set NDVI thr = 0.1 having observed that lower NDVI usually corresponded to partially snow-covered sites and or to senescent canopies (Fig. 2).The main findings of this study did not change when I varied NDVI thr in the range 0.05-0.15.As a simpler alternative to GPP int , i.e., not accounting for incoming solar radiation, I also calculated the time-integrated value of NDVI, hereafter referred to as NDVI int following: The periods from the beginning of the snow-free period to TNDVI max , hereafter referred to as P g , and from TNDVI max to the end of the first snowfall, hereafter referred to as P s , were used to decompose productivity into two components: NDVI int,g and GPP int,g , and NVI int,s and GPP int,s (Fig. 2).Note that the suffix letters g and s are used to refer to the first and the second part of the growing season, respectively.
The whole analysis was also conducted with the Enhanced Vegetation Index (Huete et al., 2002) instead of NDVI.The rationale for this alternative was to select a vegetation index which was more related to the green biomass and thus may better approximate GPP, especially during the senescence period.I did not find any significant change in the main results when using EVI.In particular, the period spanning from peak standing biomass to the first snowfall accounted for two-thirds of EVIint as is the case for NDVI int (Fig. S3a) and inter-annual variations in EVIint were of the same order of magnitude as those for NDVI int (Fig. S3b).Because results from the path analysis (see below) were also very similar with EVI-based proxies of productivity, I chose to present NDVI-based results only.

Path analysis
Path analysis represents an appropriate statistical framework to model multivariate causal relationships among observed variables (Grace et al., 2010).Here, I examined different causal hypotheses of the cascading effects of meteorological forcing, snow cover duration and phenological parameters (TNDVI max , P g , and P s ) on NDVI int and GPP int .To better contrast the processes involved during different stages of the growing season, separate models were implemented for the period of growth and the period of senescence.The set of causal assumptions is represented using directed acyclic graphs in which arrows indicate which variables are influencing (and are influenced by) other variables.These graphs may include both direct and indirect effects.An indirect effect of X1 on Y means that the effect of X1 is mediated by another variable (for example X1 → X2 → Y ).Path analysis tests the degree to which patterns of variance and covariance in the data are consistent with hypothesized causal links.To develop this analysis, three main assumptions have been made: (i) that the graphs do not include feedbacks (for example, X1 → X2 → Y → X2); (ii) that the relationships among variables can be described by linear models and (iii) that annual observations are independent, i.e., the growth response in year n is not influenced by previous years because of carryover effects.
Since I chose to focus on the inter-annual variability of growth response, I removed between-site variability by calculating standardized anomalies for each polygon.Standardized anomalies were calculated by dividing annual anomalies by the standard deviation of the time series making the magnitude of the anomalies comparable among sites.
For each causal diagram, partial regression coefficients were estimated for the whole data set and for each polygon.These coefficients measure the extent of an effect of one variable on another while controlling for other variables.Model estimates were based on maximum likelihood, and Akaike information criterion (AIC) was used to compare performance among competing models.Only ecologically meaningful relationships were tested.The model with the lowest AIC was retained as being the most consistent with observed data.
I used the R software environment (R Development Core Team, 2010) to perform all statistical analyses.Path coefficients and model fit were estimated using the package "Lavaan" (Rosseel, 2012).

Results
One hundred and twenty polygons fulfilling the selection criteria were included in the analyses.These polygons spanned 2 • of latitude and more than 1 • of longitude and were distributed across 17 massifs of the French Alps from the northern part of Mercantour to the Mont-Blanc massif (Fig. 1).Their mean elevation ranged from 1998 m to 2592 m with a median of 2250 m.Noticeably, many polygons were located in the southern and in the innermost part of the French Alps where high elevation landscapes with grassland-covered gentle slopes are more frequent, essentially because of the occurrence of flysch, a bedrock on which deep soil formation is facilitated.
A typical yearly course of NDVI and NDSI is shown in Fig. 2. The date at which the NDSI / NDVI ratio crosses the threshold of 1 was very close to the date at which NDVI crosses the threshold of 0.1.On average, NDVI max was reached 50 days after snowmelt, a period corresponding to only one-third of the length of the snow-free period (Fig. 3a).Similarly, NDVI g accounted for one-third of the NDVI int (Fig. 3b).The contribution of the first part of season was slightly higher for GPP int , though it largely remained under 50 % (Fig. 3c).Thus, the maintained vegetation greenness from TNDVI max to TSNOW fall explained the dominant contribution of the second part of the growing season to NDVIderived proxies of grassland productivity.Most of the variance in NDVI int and GPP int was accounted for by between-polygon variations that were higher during the period of senescence compared to the period of growth (Table 1).Inter-annual variations in NDVI int and GPP int represented 25 % of the total variance and were particularly pronounced at the end of the examined period, with the best year (2011) sandwiched by 2 (2010, 2012) of the 3 worst years (Fig. 4a).The two likely proximal causes of these inter- annual variations, i.e., P sf and NDVI max , showed highly contrasted variance partitioning.Between-year variation in P sf was 4 to 5 times higher than that of NDVI max (Table 1).The standardized inter-annual anomalies of P sf showed remarkable similarities with those of NDVI int and GPP int , both in terms of magnitude and direction (Fig. 4b).By contrast, the small inter-annual variations in NDVI max did not relate to inter-annual variations in NDVI int or GPP int (Fig. 4c).For example, the year 2010 had the strongest negative anomaly for both P sf and NDVI int , whereas the NDVI max anomaly was positive.There were some discrepancies between the two proxies of primary productivity.For example, the heatwave of 2003, which yielded the highest NDVI max , exhibited a much stronger positive anomaly for GPP int than for NDVI int and this was due to the unusually high frequency of clear sky during this particular summer.
The path analysis confirmed that the positive effect of the length of the period available for plant activity largely surpassed that of NDVI max to explain inter-annual variations in NDVI int and GPP int .This held true for NDVI int,g or GPP int,g -with an over-dominating effect of P g (Fig. 5a, c) -and for NDVI int,s or GPP int,s -with an over-dominating effect of P s (Fig. 5b, d).There was some support for an indirect effect of P g on productivity mediated by NDVI max , as removing the path P g → NDVI max in the model decreased its performance (Table 2).In addition to shortening the time available for growth and reducing primary productivity, a delayed snowmelt also significantly decreased the number of frost events and this had a weak positive effect on both NDVI int,g and GPP int,g (Fig. 5a, c).However, this positive and indirect effect of TSNOW melt on productivity, which amounts to (−0.46) × (−0.08) = 0.04 for NDVI int,g and (−0.46) × (−0.13) = 0.06 for GPP int,g , was small compared to the negative effect of TSNOW melt on NDVI int,g (−1 × 0.96 for NDVI int,g and −1 × 0.95 for GPP int,g ).Apart from its effect on frost events and P s , TSNOW melt also had a significant positive effect on TNDVI max with a path coefficient of 0.57, signifying that grasslands partially recover the time lost because of a long winter to reach peak standing biomass.On average, a 1-day delay in the snowmelt date translates to a 0.5-day delay in TNDVI max (Fig. S4a).
Compared to snow cover dynamics, weather conditions during the growing period had relatively small effects on both NDVI max and productivity (Fig. 5).For example, removing the effects of temperature on NDVI max and precipitation on NDVI int,g did not change model fit (Table 2).The most significant positive effects of weather conditions were observed during the senescence period and more specifically for GPP int,s with a strong positive effect of temperature (Fig. 5d).The impact of warm and dry days on incoming radiation explained why more pronounced effects of temperature and precipitation are observed for GPP int (Fig. 5d), which is dependent upon PAR (see Eq. 3), than for NDVI int (Fig. 5b).
Meteorological variables governing snow cover dynamics had a strong impact on primary productivity (Fig. 5).A warm spring advancing snowmelt translated into a significant positive effect on NDVI int,g and GPP int,g -an indirect effect which amounts to (−0.62) × (−1) × 0.95 = 0.59 (Fig. 5a, c).Heavy precipitation and low temperature in October-November caused early snowfall and shortened P s , which severely reduced NDVI int,s and GPP int,s (Fig. 5b, d).Overall, given that the senescence period accounted for twothirds of the annual productivity (Fig. 3b, c), the determinants of the first snowfall were of paramount importance for explaining inter-annual variations in NDVI int and GPP int .
Path coefficients estimated for each polygon showed that the magnitude and direction of the direct and indirect effects were highly conserved across the polygons.The climatology of each polygon was estimated by averaging growing season temperature and precipitation across the 13 years.Whatever the path coefficient, neither of these two variables explained more than 8 % of variance of the between-polygon variation (Table 3).The two observed trends were (i) a greater positive effect of NDVI max on NDVI int,g in polygons receiving more rainfall, which was consistent with the significant effect of precipitation on NDVI max (Fig. 5a) and (ii) a smaller effect of temperature and P s on GPP int,s and NDVI int,s , respectively, suggesting that the coldest polygons were less responsive to increased temperatures or lengthening of the growing period (see discussion).

Discussion
Using a remote sensing approach, I showed that inter-annual variability in NDVI-derived proxies of productivity in alpine   2).A solid line (or dotted lines) indicates a significant positive (or negative) effect at P < 0.05.Double-lined arrows correspond to fixed parameters.Abbreviations include TEMP, averaged daily mean temperature (or senescence period); PREC, averaged daily sum of precipitation; and FrEv, number of frost events.Letter g (or s) represents the initial growth period (or the senescence period), spring the months of March and April, and fall the months of October and November.grasslands was primarily governed by variations in the length of the snow-free period.As a consequence, meteorological variables controlling snow cover dynamics are of paramount importance to understand how grassland growth adjusts to changing conditions.This was especially true for the determinants of the first snowfall, given that the period spanning from the peak standing biomass onwards accounted for two-thirds of annual grassland productivity.By contrast, NDVI max -taken as an indicator of growth responsiveness -showed small inter-annual variation and weak sensitivity to summer temperature and precipitation.Overall, these results highlighted the ability of grasslands to track interannual variability in the timing of the favorable season by maintaining green tissues during the whole snow-free period and their relative inability to modify the magnitude of the growth response to the prevailing meteorological conditions during the summer.I discuss these main findings below in light of our current understanding of extrinsic and intrinsic factors controlling alpine grassland phenology and growth.
In spring, the sharp decrease of NDSI and the initial increase of NDVI were simultaneous events (Fig. 2).Previous reports have shown that NDVI may increase independently of greenness during the snow melting period (Dye and Tucker, 2003) and this has led to the search for vegetation indices other than NDVI to precisely estimate the onset of greenness in snow-covered ecosystems (Delbart et al., 2006).Here I did not consider that the period of plant activity started with the initial increase of NDVI.Instead I combined NDVI and NDSI indices to estimate the date of snowmelt and then used a threshold value of NDVI = 0.1 before integrating NDVI over time.By doing this, I strongly reduced the confounding effect of snowmelt on the estimate of the onset of greenness.That said, a remote sensing phenology may fail to accurately capture the onset of greenness for many other reasons, including smoothing procedures applied to NDVI time series, inadequate thresholds, geolocation uncertainties, mountain terrain complexity, and vegetation heterogeneity (Cleland et al., 2007;Tan et al., 2006;Dunn and de Beurs, 2011;Doktor et al., 2009).Assessing the magnitude of this error is difficult as there have been very few studies comparing ground-based phenological measurements with remote sensing data, and furthermore, most of the available studies have focused on deciduous forests (Hmimina et al., 2013;Busetto et al., 2010;but see Fontana et al., 2008).Groundbased observations collected at one high elevation site and corresponding to a single MOD09A1 pixel provide preliminary evidence that the NDVI / NDSI criterion adequately captures snow cover dynamics (Fig. S3).Further studies are needed to evaluate the performance of this metric at a regional scale.For example, the analysis of high-resolution remote sensing data with sufficient temporal coverage is a promising way to monitor snow cover dynamics in complex alpine terrain and to assess its impact on the growth of alpine grasslands (Carlson et al., 2015).Such an analysis has yet to be done at a regional scale.Despite these limitations, I am confident that the MODIS-derived phenology is appropriate for addressing inter-annual variations in NDVI int because: (i) the start of the season shows low NDVI values, and thus uncertainty in the green-up date will marginally affect integrated values of NDVI and GPP, and (ii) beyond errors in estimating absolute dates, remote sensing has been shown to adequately capture the inter-annual patterns of phenology for a given area (Fisher and Mustard, 2007;Studer et al., 2007), and this is precisely what is undertaken here.
Regardless of the length of the winter, there was no significant time lag between snow disappearance and leaf greening at the polygon level.This is in agreement with many field observations showing that initial growth of mountain plants is tightly coupled to snowmelt timing (Körner, 1999).This plasticity in the timing of the initial growth response, which is enabled by tissue preformation, is interpreted as an adaptation to cope with the limited period of growth in seasonally snow-covered ecosystems (Galen and Stanton, 1991).Early disappearance of snow is controlled by spring temperature, and our results, showing that a warm spring leads to a prolonged period of plant activity, are consistent with those originally reported from high latitudes (Myneni et al., 1997).Other studies have also shown that the onset of greenness in the Alps corresponds closely with year-to-year variations in the date of snowmelt (Stockli and Vidale, 2004) and that spring mean temperature is a good predictor of meltout (Rammig et al., 2010).This study improves upon previous works (i) by carefully selecting targeted areas to avoid mixing different vegetation types when examining growth response, (ii) by using a meteorological forcing that is more appropriate to capture topographical and regional effects compared to global meteorological gridded data (Frei and Schär, 1998), and (iii) by implementing a statistical approach en-abling the identification of direct and indirect effects of snow on productivity.
Even if there were large between-year differences in P g , the magnitude of year-to-year variations in NDVI max were small compared to that of NDVI int or GPP int (Table 1 and Fig. 4).Indeed, initial growth rates buffer the impact of interannual variations in snowmelt dates, as has already been observed in a long-term study monitoring 17 alpine sites in Switzerland (Jonas et al., 2008).Essentially, the two contrasting scenarios for the initial period of growth observed in this study were either a fast growth rate during a shortened growing period in the case of a delayed snowmelt, or a lower growth rate over a prolonged period following a warm spring.These two dynamics resulted in nearly similar values of NDVI max as TSNOW melt explained only 4 % of the variance in NDVI max (Fig. S4b).I do not think that the low variability in the response of NDVI max to forcing variables is due to a limitation of the remote sensing approach.First, there was a high between-site variability of NDVI max , indicating that the retrieved values were able to capture variability in the peak standing aboveground biomass (Table 1).Second, the mean NDVI max of the targeted areas is around 0.7 (Fig. 4b), i.e., in a range of values where NDVI continues to respond linearly to increasing green biomass and Leaf Area Index (Hmimina et al., 2013).Indeed, many studies have shown that the maximum amount of biomass produced by arctic and alpine species or meadows did not benefit from the experimental lengthening of the favorable period of growth (Kudo et al., 1999;Baptist et al., 2010), or to increasing CO 2 concentrations (Körner et al., 1997).Altogether, these results strongly suggest that intrinsic growth constraints limit the ability of high elevation grasslands to enhance their growth under ameliorated atmospheric conditions.More detailed studies will help us understanding the phenological response of different plant life forms (e.g., forbs and graminoids) to early and late snow-melting years and their contribution to ecosystem phenology (Julitta et al., 2014).Other severely limiting factors -including nutrient availability in the soil -may explain this low responsiveness (Körner, 1989).For example, Vittoz et al. (2009) emphasized that year-to-year changes in the productivity of mountain grasslands were primarily caused by disturbance and land use changes that affect nutrient cycling.Alternatively, one cannot rule out the possibility that other bioclimatic variables could better explain the observed variance in NDVI max .For example, the inter-annual variations in precipitation had a slight though significant effect on NDVI max (Fig. 5a, c) suggesting that including a soil-water balance model might improve our understanding of growth responsiveness, as suggested by Berdanier and Klein (2011).
Many observations and experimental studies have also pointed out that soil temperature impacts the distribution of plant and soil microbial communities (Zinger et al., 2009), ecosystem functioning (Baptist and Choler, 2008), and flowering phenology (Dunne et al., 2003) et al., 2008;Freppaz et al., 2007).Thus, an improvement of this study would be to test, not only for the effect of presence/absence of snow, but also for the effect of snowpack height and soil temperature on NDVI max and growth responses of alpine pastures.Regional climate downscaling of soil temperature at different depths is currently under development within the SAFRAN-CROCUS-MEPRA model chain and there will be future opportunities to examine these linkages.Nevertheless, the results showed that at the first order the summer meteorological forcing was instrumental in controlling GPP int,s , without having a direct effect on NDVI max (Fig. 5b, d).In particular, positive temperature anomalies and associated clear skies had significant effects on GPP int,s .Moreover, path analysis conducted at the polygon level also provided some evidence that responsiveness to ameliorated weather conditions was less pronounced in the coldest polygons (Table 3), suggesting stronger intrinsic growth constraints in the harshest conditions.Collectively, these results indicated that the mechanism by which increased summer temperature may enhance grassland productivity was through the persistence of green tissues over the whole season rather than through increasing peak standing biomass.The contribution of the second part of the summer to annual productivity has been overlooked in many studies (e.g., Walker et al., 1994;Rammig et al., 2010;Jonas et al., 2008;Jolly et al., 2005) that have primarily focused on early growth, or on the amount of aboveground biomass at peak productivity.Here, I showed that the length of the senescing phase is a major determinant of inter-annual variation in growing season length and productivity and hence that temperature and precipitation in October-November are strong drivers of these inter-annual changes (Fig. 5b, d).The importance of autumn phenology was recently re-evaluated in remote sensing studies conducted at global scales (Jeong et al., 2011;Garonna et al., 2014).A significant long-term trend towards a delayed end of the growing season was noticed for Europe and specifically for the Alps.In the European Alps, temperature and moisture regimes are possibly under the influence of the North Atlantic Oscillation (NAO) phase anomalies (Beniston and Jungo, 2002) in late autumn and early winter.This opens the way for research on teleconnections between oceanic and atmospheric conditions and the regional drivers of alpine grassland phenology and growth.
Eddy covariance data also provided direct evidence that the second half of the growing season is a significant contributor to the annual GPP of mountain grasslands (Chen et al., 2009;Rossini et al., 2012;Li et al., 2007;Kato et al., 2006).However it has also been shown that while the combination of NDVI and PAR successfully captured daily GPP dynamics in the first part of the season, NDVI tended to provide an overestimate of GPP in the second part (Chen et al., 2009;Li et al., 2007).Possible causes include decreasing light-use ef-ficiency in the end of the growing season in relation to the accumulation of senescent material and/or the "dilution" of leaf nitrogen content by fixed carbon.Noticeably, the main findings of this study did not change when NDVI was replaced by EVI, a vegetation index which is more sensitive to green biomass and thus may better capture primary productivity.Consistent with this result, Rossini et al. (2012) did not find any evidence that EVI-based proxies performed better than NDVI-based proxies to estimate the GPP of a subalpine pasture.Further comparison with other vegetation indexes -like the MTCI derived from MERIS products (Harris and Dash, 2010) -will contribute to better evaluations of NDVI-based proxies of GPP.
A strong assumption of this study was to consider that the LUE parameter is constant across space and time.There is still a vivid debate on the relevance of using vegetation specific LUE in remote sensing studies of productivity (Yuan et al., 2014;Chen et al., 2009).Following Yuan et al. (2014) I have assumed that variations in light-use efficiency are primarily captured by variations in NDVI because this vegetation index correlates with structural and physiological properties of canopies (e.g., leaf area index, chlorophyll, and nitrogen content).Multiple sources of uncertainty affect remotely sensed estimates of productivity and it is questionable whether the product NDVI times PAR is a robust predictor of GPP in alpine pastures.The estimate of absolute values of GPP and its comparison across sites was not the aim of this study that focuses on year-to-year relative changes of productivity for a given site.It is assumed that limitations of a light-use efficiency model are consistent across time and that these limitations did not prevent the analysis of the multiple drivers affecting inter-annual variations in remotely sensed proxies of GPP.At present, there is no alternative for regional-scale assessment of productivity using remote sensing data.In the future, possible improvements could be made by using air-borne hyperspectral data to derive spatial and temporal changes in the functional properties of canopies (Ustin et al., 2004) and assess their impact on annual primary productivity.

Conclusions
I have shown that the length of the snow-free period is the primary determinant of remote sensing-based proxies of primary productivity in temperate mountain grasslands.From a methodological point of view, this study demonstrated the relevance of path analysis as a means to decipher the cascading effects and relative contributions of multiple predictors on grassland phenology and growth.Overall, these findings call for a better linkage between phenomenological models of mountain grassland phenology and growth and land surface models of snow dynamics.They open the way to a process-based, biophysical modeling of alpine pastures growth in response to environmental forcing, follow-ing an approach used in a different climate (Choler et al., 2010).Year-to-year variability in snow cover in the Alps is high (Beniston et al., 2003) and climate-driven changes in snow cover are on-going (Hantel et al., 2000;Keller et al., 2005;Beniston, 1997).Understanding the factors controlling the timing and amount of biomass produced in mountain pastures thus represents a major challenge for agro-pastoral economies.
The Supplement related to this article is available online at doi:10.5194/bg-12-3885-2015-supplement.
Figure 1.(a) Location map of the 121 polygons across the 17 climatologically defined massifs of the French Alps.(b) Number of polygons per massif.

Figure 2 .
Figure 2.Yearly course of NDVI and NDSI showing the different variables used in this study: date of snowmelt (TSNOW melt ), maximum NDVI (NDVI max ) and date of NDVI max (TNDVI max ), date of snowfall (TSNOW fall ), length of the snow-free period (P sf ), length of the initial growth period (P g ), length of the senescence period (P s ), and time-integrated NDVI over the growth period (NDVI int,g ) and over the senescence period (NDVI int,s ).

Figure 3 .
Figure 3. Frequency distribution of the relative contribution of P g and P s to P sf (a), of NDVI int,g and NDVI int,s to NDVI int (b), and of GPP int,g and GPP int,s to GPP int (c).Values were calculated for each year and for each polygon.

FrEvFigure 5 .
Figure 5. Path analysis diagram showing the interacting effects of meteorological forcing, snow cover duration, and NDVI max on NDVI int (a, b) and GPP int (c, d).For each proxy of productivity, separate models for the period of growth (a, c) and the period of senescence (b, d) are shown.Line thickness of arrows is proportional to standardized path coefficients which are indicated on the right or above each arrow.Values in italics indicate paths that can be removed without penalizing model AIC (see Table2).A solid line (or dotted lines) indicates a significant positive (or negative) effect at P < 0.05.Double-lined arrows correspond to fixed parameters.Abbreviations include TEMP, averaged daily mean temperature (or senescence period); PREC, averaged daily sum of precipitation; and FrEv, number of frost events.Letter g (or s) represents the initial growth period (or the senescence period), spring the months of March and April, and fall the months of October and November.

Table 1 .
Variance partitioning into between-polygon and between-year components for the set of predictors and growth responses included in the path analysis.

Table 2 .
Model fit of competing path models.AIC is the Akaike information criteria value and AIC is the difference in AIC between the best model and alternative models.

Table 3 .
Relationships between mean temperature or precipitation of polygons and the path coefficients estimated at the polygon level.Only significant relationships are shown.

www.biogeosciences.net/12/3885/2015/ Biogeosciences, 12, 3885-3897, 2015 P. Choler: Growth response of grasslands to snow cover duration ing
. More specifically, the lack of snow or the presence of a shallow snowpack dur-winter increases the frequency of freezing and thawing events with consequences on soil nutrient cycling and aboveground productivity (Kreyling