Frédéric Marin
email: frederic.marin@utc.fr

Kevin Lepetit

Laetitia Fradet

Clint Hansen

Ben Khalil

Mansour

Using accelerations of single inertial measurement units to determine the intensity level of light-moderate-vigorous physical activities: Technical and mathematical considerations

Keywords: Physical activity, gait, run, IMU, Wearable inertial sensor, motion capture, biomechanics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Physical activity improves the general well-being [START_REF] Kubota | Physical Activity and Lifetime Risk of Cardiovascular Disease and Cancer[END_REF]. However, for an efficient physical activity program, the optimal physical activity intensity has to be determined. Planning adequate physical activity and then prescribing an optimal intensity for those is a challenge for health care policies (WHO, 2010). The relevance of on-body sensors to monitor physical activity is wellestablished [START_REF] Marschollek | Associations between sensor-based physical activity behaviour features and healthrelated parameters[END_REF]) and inertial measurement units (IMUs) have become a popular technical solution [START_REF] Rault | A survey of energy-efficient context recognition systems using wearable sensors for healthcare applications[END_REF] within their constraints [START_REF] Kerr | Comparison of Accelerometry Methods for Estimating Physical Activity[END_REF].

Although multiple studies have dealt with physical activity monitoring [START_REF] Fradet | Classification of physical activities based on body-segments coordination[END_REF], there has been little investigation of the mathematical and technical rationale regarding the methodology used to establish the correspondence between the inertial sensors signal (commonly acceleration) and the physical activity measure itself. The technical rationale or the mathematical background are most of the time superficial, not demonstrated and, in most cases, only a final output without any technical background is provided. In addition, as far as the technical implementation is concerned, there is no consensus as to sensor placement or input data prerequisites [START_REF] Nez | Identification of Noise Covariance Matrices to Improve Orientation Estimation by Kalman Filter[END_REF]. Often various positions are suggested for the use of the sensor [START_REF] Sandroff | Accelerometer measured physical activity and the integrity of the anterior visual pathway in multiple sclerosis[END_REF] even though the hip is already recognized as providing better results than the wrist (Cleland et al. 2013). Hardware heterogeneity [START_REF] Sandroff | Accelerometer measured physical activity and the integrity of the anterior visual pathway in multiple sclerosis[END_REF] is another well-known issue, potentially causing different results [START_REF] Nez | Comparison of calibration methods for accelerometers used in human motion analysis[END_REF]). An example is the varying sampling rate according to the hardware of the sensor (Yang 2010) and the difference between continuous and intermittent data collection. Some recent studies have demonstrated the influence of sampling frequency on PA estimation [START_REF] Brønd | Sampling frequency affects the processing of Actigraph raw acceleration datato activity counts[END_REF].

Regarding the mathematical background, the most common outputs are "activity count", identified over a specific period of time [START_REF] Yang | A review of accelerometry-based wearable motion detectors for physical activity monitoring[END_REF]. However, the exact calculation of such parameters is often missing, and leaves room for speculation. Most of the time, this activity count derives from acceleration measurements of the IMU and analytical approaches have been proposed (Chen et al. 2005). The acceleration signal is rectified and integrated by a user-defined epoch. An "activity" is counted once the resulting signal exceeds a certain threshold. However, this threshold is often undocumented.

To quantify movement performance various metrics based on IMU have been proposed in the literature [START_REF] Lepetit | Optimized scoring tool to quantify the functional performance during the sit-to-stand transition with a magneto-inertial measurement unit[END_REF]. Some are easy to determine, such as the vector magnitude minus one (Van Hees 2013 ), others, however, require numerous calculi [START_REF] Fradet | Classification of physical activities based on body-segments coordination[END_REF]. After this metric computation, regression equations are proposed to estimate energy expenditure [START_REF] Rothney | Validation of the ActiGraph tworegression model for predicting energy expenditure[END_REF] or cut-points are proposed for sedentary time and physical activity intensity classification [START_REF] Migueles | Accelerometer Data Collection and Processing Criteria to Assess Physical Activity and Other Outcomes: A Systematic Review and Practical Considerations[END_REF]. The level of activity as "light", "moderate" or "vigorous" is based on Metabolic Equivalent of Task (MET) [START_REF] Ainsworth | Compendium of Physical Activities: a second update of codes and MET values[END_REF] even though the relationship between IMU based metrics and MET is questionable [START_REF] Yang | A review of accelerometry-based wearable motion detectors for physical activity monitoring[END_REF].

Previous elements showed that the use of a single IMU to qualify and classify the intensity of physical activity did not find consensus in terms of set up and mathematical rationale for IMU deduced metrics. Consequently, IMU's users are faced with multiple options to place the IMU on the body as well as with various calculations to score the level of activity. In addition, in some propriety systems, this calculation is hidden. In this context, the aim of this paper is (i) to provide direct comparison of the efficiency of various IMU deduced metrics according to IMU localization to quantify and differentiate light-moderate-vigorous physical activities, and (ii) to assess the relationship between those metrics and the MET.

Material and Methods

147 participants (108 females and 39 males; age: 21.0 years (SD 2.0 y.), body mass index: 21.7 kgm -2 (SD 2.6 kgm -2 ) voluntarily participated in the experiment after signing a statement of informed consent pertaining to the experimental procedure as required by the Helsinki declaration. Data were collected (figure 1) in two sessions: one when the IMU (Opal, APDM) was located at the wrist and the second one when the IMU was located at the ankle. All sensors were calibrated before each session to avoid sensor drift or offsets according to the manufacturer's procedure. The sample rate was 128 Hz. Participants were asked to maintain a standing posture for 2 minutes, followed by self-pace walking (466 m), followed by another 2 minutes standing, followed by a self-pace run over (700 m).

The session ended in a quiet standing position.

According to the physical activity classification (WHO 2010), the standing position was defined as light physical activity (L), walking as moderate physical activity (M) and running as a vigorous one (V).

In order to individually quantify the level of the physical activity, the corrected MET was computed according to the physical activity, age, gender, height, weight and the locomotion speed [START_REF] Ainsworth | Compendium of Physical Activities: a second update of codes and MET values[END_REF].

The post processing of the IMU's measurements was computed in five steps and this study focuses specifically on the accelerometer data (Figure 2):

a) Accelerations deduced from IMU's measurements

According to the hardware of the sensors, IMU measurements are finite data with a frequency of herz . Consequently, if the session time is T in seconds, we obtained data with n measurements (eq. 1).

= .

(eq. 1)

In a second step, three different accelerations were computed: raw accelerations, IMU accelerations and body accelerations.

The raw sensor acceleration data is expressed in the reference frame of the sensor, named at each measure ∈ 1, , it corresponds to a vector in ℝ according to the values on the 3 axis of the sensors in the case of 3D accelerometers (eq. 2).

= (eq. 2)

As the sensor is an inertial measurement unit including a 3D gyroscopes and a 3D magnetometers it is possible to determine the orientation matrix

!"#$
For each measure describing the orientation of the IMU relative to the laboratory frame using data fusion algorithms [START_REF] Nez | Identification of Noise Covariance Matrices to Improve Orientation Estimation by Kalman Filter[END_REF].

The accelerations were then expressed in the laboratory reference frame %&' (eq. 3). Classically, the

x axis and y axis of the laboratory are in the horizontal plane and the z axis is vertical (eq. 4).

%&'

= %&' (eq. 3)

%&' = %&' %&' %&'
(eq. 4)

Due to the microelectromechanical components, the measure of the acceleration includes the acceleration of the rigid body (i.e. the IMU) and the gravity [START_REF] Tamura | Chapter 2.2 Wearable Inertial Sensors and Their Applications in book Wearable Sensors : Fundamentals, Implementation and Applications[END_REF]. Accordingly the z axis of the laboratory reference frame is vertical such as the gravitational component ( (1 ( = 9.806 ms -2 ) can be subtracted ( eq. 5 ). We deduced then the acceleration of the rigid body in the laboratory reference frame, named $)* ( eq. 6 )

$)* = %&' -( = %&' %&' %&' -( (eq. 5) $)* = $)* $)* $)* (eq. 6) b)

Physical activities metrics calculation

To quantify the intensity of physical activities at each i-th measurement, various parameters deduced from the accelerations previously obtained are proposed in the literature.

Van Hees (2013) and [START_REF] Hildebrand | Age group comparability of raw accelerometer output from wrist-and hip-worn monitors[END_REF] proposed the "vector magnitude minus one" named ENMO i (eq. 7).

1234 = 5 6 + 6 + 6 -( (eq. 7)

Vähä-Ypyä (2015) proposed to compute the "amplitude deviation" named AD (eq 8).

:

; = <5 6 + 6 + 6 - 1 = 5 > 6 + > 6 + > 6 " ?@A < (eq. 8)
Based on $)*

, Taruma (2014) suggested the use of the" vertical acceleration", named VA (eq. 9).

C: = $)*
(eq. 9)

VanHees (2011) identified the "vector magnitude" of the acceleration named VM (eq.10) as an indicator of the intensity of the physical activity.

C3 = D E$)* 6 + E$)* 6 + E$)* 6
(eq.10)

Chen (1997) suggested a direct estimation of the energy expenditure named EEact ) (eq. 11) based on the magnitude of the horizontal acceleration and vertical acceleration. For each subject (S), and for the wrist (W) and ankle session (A), sensor data were segmented to extract the metric MPA from L, M and V physical activity [START_REF] Fradet | Classification of physical activities based on body-segments coordination[END_REF]. At this stage, MPA is a ℝ " vectors. In addition, is subject dependent because duration walking and running sessions were selfpaced.

11 IJ = K LD E$)* 6 + E$)* 6 M NA + O K P E$)* Q N6 with, R1 = 2.
Hammerla (2013) suggest analyzing the quantile function of the metric M PA. The quantile function a I is defined as the inverse of the empirical cumulative distribution. a I is a strict monotonically increasing function which associates a cumulative probability c from zero to one with the level of the value of the metric MPA [START_REF] Gilchrist | Interpreting regression models in clinical outcome studies[END_REF]. Consequently, c = 0 would be associated with min MPA , and c = 1 with max MPA (figure 4). Then, a I was linearly interpolated (0 to 1 with a step of 0.1). We obtained discrete quantile function d & ef as ℝ AgA vector for a metrics MPA associated at a subject, a sensor localization and session.

c) Data reductions

Data reduction consist in reducing data dimension. It consisted in reducing MPA which is ℝ " vectors to ℝ h vector with i < .

Usually in literature, data reduction of a metric 3 kl ∈ ℝ " to a single value in ℝ is concerned, authors suggest to determine the mean value (MEAN) (Vähä-Ypyä 2015) and the root mean square (RMS) value (Taruma 2014).

For the data reduction of a quantile function, L-moments method as higher order descriptive statistical tool is used (Karvanen 2005) and computed as follows: First, λ-coefficients and L-moments are computed with numerical integration performed by rectangle rule.

m s = n 30I 6 -30I 6 + 12I -1 a I oI A g ≈ = q30 L > -1 101 M -30 L > -1 101 M 6 + 12 L > -1 101 M -1r d & ef
AgA ?@A (eq. 17)

L-moments analogous to the coefficient of variation (CV), the skewness of the distribution and the coefficient of kurtosis could be deduced (Karvanen 2005):

t uv = m 6 m A (eq .18) t #w! "!## = m m 6
(eq .19) To quantify the similarity of the inter-subject of the quantile curves, the coefficient of multi correlation (CMC) on the quantile functions d & ef of all parameters for all subjects and for each physical activity and sensor's positions was computed in order to analyze the inter-subject reliability of this analysis (Kabada 1989). CMC scores the similarity between curves by a value between 0 (very different) to 1 (very similar) (Kabada 1989).

t
Finally, the relationship between data reductions MEAN and RMS and MET values of the subject was investigated by a linear regression model with a robust linear least mean square method. The relevance of the linear regression model by the coefficient of determination R 2 was reported [START_REF] Draper | Applied regression analysis[END_REF]. All computation and statistics analyses were performed in Matlab (R2016a, The Mathworks Inc, Natick, MA).

Results

Based on the Kruskal-Wallis test, we noticed no significant difference between the wrist and ankle sensor for the mean values of the ENMO, the RMS values of the ENMO and AD metrics for only for the light physical activities (Figures 3). Consequently, for all other parameters, we noticed a significant difference according to the localization of the sensor and the level of the physical activities (Figure 3). When the sensor is located on the wrist, we noticed that the box plots for all data reductions during light and moderate activities were close, which means that ranges of data reductions for light activity were close to that during a moderate one. These observations presented a lower differentiation ability. On the contrary, box plots for data reductions during vigorous activity were more clearly separated from those occurring during light and moderate ones, which suggests that vigorous level of physical activity could easily be differentiated.

For the quantile function of all metrics (Figure 4), we found a CMC value from 0.74 to 0.98 which demonstratred as a good to very good inter-subject reliability [START_REF] Rezgui | On the imitation of CP gait patterns by healthy subjects[END_REF]. The values of less than 0.80 were obtained for a light level of physical activity at the wrist localization for the AD, VM and EEact metrics and at the ankle localization for the VA and EEact metrics. High level of CMC ( > 0.9) for the three levels of activity were obtained for the quantile function of the metrics ENMO when the sensor is located at the wrist or ankle, and only for the ankle localization for the quantile function of metrics AD, VA and VM at moderate and vigorous physical activities.

Figure 5 showed results for L-CV, L-Skewness and L-Kurtosis of the quantile function of all metrics.

We noticed no significant difference between the wrist and ankle sensor for the L-CV of VA and EEact at light activity,and for L-Skewness of VM and L-Kurtosis of ENMO and EEact at vigorous activity .

We also found no significant difference between the light and the moderate activities for the L-CV of ENMO at the ankle, and between the moderate and vigorous activities for the L-CV of AD and VM at the wrist, for the L-CV of AD at the ankle, and for the L-Skewness of AD and EEact of the ankle. We also noticed the box plots for L-Skewness and L-Kurtosis of light level are more clearly separated from moderate and vigorous ones, which suggests that light level of physical activity could easily be differentiated from the moderate and vigorous ones using these metrics.

Good to very good linear relationship between the MEAN and RMS of all parameters with the MET values has been shown with values of the coefficient of determination R 2 from 0.78 to 0.97 (Figure 6). Indeed, based on the wrist sensor data, 97% of the linear relationships between the mean values and AUC of the AD and VM with the MET values has been accounted for. For the ankle sensor data, 97% of the linear relationships between the mean values of the ENMO with the MET. We also noticed visible difference on the slope obtained for wrist and ankle based data for MEAN-Met linear function of the metrics ENMO and VM (figure 6).

Discussion

The purpose of this study was to highlight the influence of IMU deduced metrics to quantify three levels of physical activity intensity and to compare between sensor locations. Our findings demonstrate that the five data reductions (MEAN, RMS, L-CV, L-Skewness and L-Kurtosis) of five metrics ENMO, AD, VA VM ,EEact) at the two localizations (wrist and ankle) could mainly differentiate between three levels of physical activity (light, moderate or vigorous). The second purpose was to test the relationship between those metrics and the MET, for which we found strong correlations and high coefficients of determination.

Research on post-processing methods of IMU based accelerometer data to quantify and classify physical activity is abundant [START_REF] Fradet | Classification of physical activities based on body-segments coordination[END_REF]) , but no standard procedure or reference data are available to test the reliability and the robustness of these methods. While some metrics like ENMO (eq. 7) and AD (eq. 8) are easy to compute (only based on 3D accelerations), obtaining metrics for VA (eq. 9) and VM (eq. 10) require additional sensors (e.g. gyroscope / magnetometer) and an additional sensor fusion algorithm [START_REF] Nez | Identification of Noise Covariance Matrices to Improve Orientation Estimation by Kalman Filter[END_REF]. Concerning the MET modeling, the metric EEACT (eq. 11) links the IMU measures and the metabolic energy using nonlinear equations. Non-linear equations present a computational cost and their generalizability to the whole population remains questionable.

Our results show a strong relationship between the MET estimation and the computed metrics (R 2 >0.8) and are consistent with previous research (Hamilton et al. 2015). This emphasizes the relevance of acceleration based metrics when estimating the level of activity [START_REF] Crouter | Estimating energy expenditure using accelerometers[END_REF] even though reliability might be low (Hargen et al. 2017). The low reliability is due to the counts.min -1 (cpm) is a common metric to quantify physical activity based on the acceleration (Aspvik.et al. 2016).

Cpm values represent the number of how many times the magnitude of acceleration reaches a defined threshold [START_REF] Migueles | Accelerometer Data Collection and Processing Criteria to Assess Physical Activity and Other Outcomes: A Systematic Review and Practical Considerations[END_REF]) and consequently does not allow to quantify the level of acceleration. This metric is often a proprietary format by the IMU manufacturers and exact calculation and processing is unknown [START_REF] Paul | Comparison of two different physical activity monitors[END_REF] leading to a heterogeneity of the cpm values [START_REF] Alhassan | Accuracy of accelerometer regression models in predicting energy expenditure and METs in children and youth[END_REF]. This emphasizes the importance of metric based computations that include the magnitude of the acceleration to quantify the intensity of physical activity (van Hees 2011). We also noticed a higher correlation with the metrics based on the wrist sensor, which could be linked to the relation between the swinging of the arms and the energy efficiency of the locomotion [START_REF] Arellano | The effects of step width and arm swing on energetic cost and lateral balance during running[END_REF] We focused our investigation on the IMU located on the wrist and ankle. The rational of this choice is the "wearability" [START_REF] Gemperle | Design for Wearability, Digest of Papers[END_REF] of the wrist and ankle, which are convenient for the subject and easy to fix an IMU with adjustable belt. Concerning the sensor location, the computed metrics showed different values and emphasized the importance of the sensor location [START_REF] Yang | A review of accelerometry-based wearable motion detectors for physical activity monitoring[END_REF] depending on the type of physical activity investigated [START_REF] Altini | Estimating energy expenditure using body-worn accelerometers: a comparison of methods, sensors number and positioning[END_REF]. The current study excluded the waist localization due to the risk of wobbling mass i.e. soft tissue artifacts [START_REF] Camomilla | A Human movement analysis: The soft tissue artefact issue[END_REF] affecting the acceleration measurement.

The data reduction techniques in this study ranged from rather descriptive parameters (MEAN, RMS) to quantile function. The quantile function presents the advantage to preserve the whole distribution of the accelerations, which is valuable for physical activity recognition [START_REF] Hammerla | On preserving statistical characteristics of accelerometry data using their empirical cumulative distribution ISWC '13[END_REF]. In order to describe this distribution, we used the L-moment formalism [START_REF] Hosking | L-moments_ analysis and estimation of distributions using linear combinations of order statistics[END_REF] which is known for its robustness. L-moments such as skewness and kurtosis are part of the high-order statistics (HOS) parameters [START_REF] Karvanen | Estimation of quantile mixtures via L-moments and trimmed L-moments[END_REF]) and allow for the description of the statistical shape distribution (Boudaoud et al. 2009). In addition, they could be used as parameters in a classification procedure [START_REF] Jones | Do changes in feature selection parameters influence the classification of knee rehabilitation exercises when using body worn accelerometer data?[END_REF]. We refrained from computing as data reduction the integral per unit time [START_REF] Bouten | Assessment of energy expenditure for physical activity using a triaxial accelerometer[END_REF][START_REF] Tsurumi | Estimation of Energy Expenditure during Sedentary Work with Upper Limb Movement[END_REF][START_REF] Sandroff | Accelerometer measured physical activity and the integrity of the anterior visual pathway in multiple sclerosis[END_REF] or the quadratic integral per unit time [START_REF] Lu | Sports biomechanics infomation acquisition and analysis based on Miniature Inertial Measurement[END_REF] as the numerical integration by rectangle rule is associated at a summation meaning that the integral per unit time is homogenous to the mean value and the quadratic integral to the RMS. Our result suggested that, on the one hand MEAN and RMS data reductions of all metrics present a better ability to largely discriminate between vigorous level of physical activity as opposed to low and moderate ones, on the other hand, L-skewness and L-kurtosis are more efficient to differentiate between low-level against moderate and vigorous activity. The combination of these observations could be a perspective to better segment low, moderate and vigorous levels during a continuous physical activity.

Our study also presented limitations. The reliability of the measurement of accelerations is a limitation as the calibration and sensor shift can have a dramatic effect on the accuracy [START_REF] Ricci | On the Orientation Error of IMU: Investigating Static and Dynamic Accuracy Targeting Human Motion[END_REF], Nez et al. 2017). However as the duration of the measurements was relatively short (less than 20min), we expected that the effect of such phenomena are limited and longer measurements will require drift correction (Nez et al. 2017). The investigation of walking and running as moderate and vigorous physical activities could be questionable, especially as numerous other moderate and vigorous physical activities do exist. However, walking and running are the most common and accessible physical activities for a large part of the population (WHO 2010) and way of monitoring physical activity is still controversial [START_REF] Aadland | Reproducibility of objectively measured physical activity: Reconsideration needed[END_REF]). In addition, our population is homogeneous in terms of age and health status, which cannot reflect all variations in walking and running matters. Our study reflected an "ideal situation" and could be seen as a baseline for further specific investigation of a more heterogeneous population, including children and senior citizens, as well as a pathological population.

Conclusion

When it comes to quantifying a physical activity, there is currently a lot of heterogeneity in terms of computation cost and output, namely as far as metrics and data reduction approaches are concerned. Our findings suggest a good method for differentiation applied to standing, walking, and running, all three popular physical activities of a larger population. According to the various metrics available to quantify and analyze physical activity and differences between them, our recommended "best practice" would be to fully document methods of calculation and set-up. Our results also suggest that MEAN and RMS data reductions are more suitable to discriminate between a vigorous level of physical activity as opposed to low and moderate ones. On the other hand, L-skewness and Lkurtosis are better at differentiating between low-level versus moderate or vigorous activities. These findings can be implemented in classification algorithms.

To conclude, wearable devices as a preview of "mobile monitoring of neuro-musculoskeletal performance" confirm a high ability for monitoring biomechanical behaviours of healthy subjects during their daily activities, but need rigorous methods and practice to develop fully their potential.

  Kolmogorov-Smirmov tests were performed on all data reductions. All null hypotheses of a standard normal distribution at the 5% significance level were rejected. It motivated to do the Kruskal-Wallis test, a non-parametric version of the classical one-way ANOVA to compare data between each other.Then, a Kruskal-Wallis test between the data reductions (MEAN, RMS, L-CV, L-Skewness and L-Kurtosis) of the metric MPA of a single localization (wrist or ankle) was performed to compare the significant difference between the level of activity. A p-values of less than 0.001 were rejected as the null hypothesis. The null hypothesis is defined when median values of the data reduction of the metric MPA for L and M levels , or, for M and V levels of physical activity are equal.

Figure 2 :Figure 3 :

 23 Figures captions

Figure 5 :

 5 Figure 5 : Descriptive statistics of the L-CV, L-Skewness and L-Kurtosis values of all metrics. The sign "≠" refers to no significant difference between the pointed values.

Figure 6 :

 6 Figure 6 : Graph of the relationship between the MEAN and RMS values of all metrics with the MET evaluation of the standing, walking and running activities.

  At this stage, five PA metrics are available; ENMO, AD, VA, VM and EEact which are ℝ " vectors were defined as indicators of the level of physical activity based on accelerations measured or deduced by IMU. The metric MPA is one of the five PA metrics: MPA ∈{ ENMO, AD, VA VM ,EEact }

	, for female subjects , and for male subjects	K = O K = O K = 38.903 #U ?!VW -672.38 12.813 #U ?!VW + 843.22 1000 1000 38.903 #U ?!VW -10.06 1000	(eq. 11)
	(3 #U ?!VW is the mass of the subject in kg)
		R2 =	663 #U ?!VW + 146.72 1000 -3.853 #U ?!VW + 968.28 1000

Conflict of interest statement

The authors have no conflicts of interest to report related to this study.


Karl Grigulis
email: karl.grigulis@univ-grenoble-alpes.fr

Sandra Lavorel

Simple field-based surveys reveal climate-related anomalies in mountain grassland production

Keywords: grassland biomass, climate response, long-term monitoring, tipping-point, mountain livestock farming

Sensitivity of grassland biomass production to climate is critical to impacts on multiple ecological processes and ecosystem services. Understanding its climate determinants is essential for climate change adaptation. This requires long-term monitoring, using robust methods that are appropriated by stakeholders. We tested the sensitivity of easily measured sward height to interannual climate variation in mountain grasslands. Using twelve consecutive years of measurements across 67 grassland plots representative of six land-use types associated with different landscape positions, we show that peak green biomass increased with mean summer months (June and July) maximum temperature. Different land-use types responded to specific combinations of climate parameters, but all except higher-elevation summer pastures were sensitive to summer months temperatures. We did not detect any effects of drought, with summer precipitation instead decreasing peak biomass of some grasslands due to cooling and cloudiness, also suggesting that soil water recharge from snowmelt was enough to sustain the first growth cycle. Summer pasture peak biomass decreased with number of frosts during the onset of growth in May. These result support the robustness and sensitivity of sward height as an indicator for climate response of peak fodder biomass. Differential responses across land-use types suggest some resource complementarity which can support tactical adaptation for farmers. During the three recent hottest summers (2015, 2017 and 2018) production was well below predicted values from actual temperatures, suggesting a potential regime shift when the vegetative growth period is shortened by temperature-driven acceleration in phenology and/or heat stress combined with high light intensity causing physiological damage. The baseline regime and the anomalies in hottest years need confirmation for longer time series and across a greater geographic extent. Further effects of drought and of an extended growing season are also likely for post-harvest or grazing regrowth.

Introduction

Grasslands are a critical asset in mountains, where they provide multiple benefits to societies and support large biological and cultural diversity [START_REF] Körner | Mountain Systems[END_REF]. In the European Alpine Space they cover 18% of the total area and supply an estimated average of 44.8 M tons of dry matter per hectare and per year, supporting livestock for dairy and meat production [START_REF] Jäger | Mapping of supply and demand of biomass from alpine grassland[END_REF][START_REF] Schirpke | Integrating supply, flow and demand to enhance the understanding of interactions among multiple ecosystem services[END_REF]. Fodder production is however highly variable both spatially and temporally [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF][START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF][START_REF] Schirpke | Integrating supply, flow and demand to enhance the understanding of interactions among multiple ecosystem services[END_REF]. In particular, effects of interannual climate variability, extreme events and trend change on fodder production will have considerable consequences on the futures of mountain agriculture and landscapes [START_REF] Darnhofer | Resilience and why it matters for farm management[END_REF][START_REF] Nettier | Resilience as a framework for analysing the adaptation of mountain summer pasture systems to climate change[END_REF][START_REF] Sautier | Assessing climatic exposure of grasslandbased livestock systems with seasonal-scale indicators[END_REF].

Biomass production from temperate grasslands is controlled for a given composition and management by incident radiation, temperature, soil moisture (determined by the balance between precipitation and evapotranspiration for a given soil) and available nutrients (Duru et al., 2010a;[START_REF] Fraser | Production changes in response to climate change[END_REF]. In cold biomes like mountains, snow cover duration is also critical for determining grassland production potential through vegetation composition (Carlson et al., 2015;[START_REF] Carlson | Observed long-term greening of alpine vegetation-a case study in the French Alps[END_REF]Liu et al., 2018a), and through interannual variability of production [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF][START_REF] Liu | Deciphering impacts of climate extremes on Tibetan grasslands in the last fifteen years[END_REF]. Studies at national to continental scales have linked inter-annual variations of grassland production to simple meteorological variables (temperature, precipitation) at varying time resolutions (e.g. monthly or sub-monthly means) (reviewed by [START_REF] Fraser | Production changes in response to climate change[END_REF]). Further composite variables such as growing degree days (the sum of positive mean daily temperatures until harvest or specific phenological states), potential evapotranspiration and its difference with precipitation, or number of frost days have also been used for modelling observed variations [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF][START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF] and for mechanistic models of production (Calanca et al., 2016;Duru et al., 2010a).

The understanding of climate responses of grassland biomass production is supported by combinations of observations, experiments and modelling [START_REF] Fraser | Production changes in response to climate change[END_REF]. Long-term monitoring offers a unique opportunity for detecting response patterns, parameterising and testing models [START_REF] Rogora | Assessment of climate change effects on mountain ecosystems through a cross-site analysis in the Alps and Apennines[END_REF][START_REF] Trnka | A simple statistical model for predicting herbage production from permanent grassland[END_REF]. Low-cost monitoring is also a means to empower stakeholders by implicating them in data collection, interpretation and communication [START_REF] Dobremez | Les alpages sentinelles : un dispositif original pour une nouvelle forme de gouvernance partagée face aux enjeux climatiques[END_REF][START_REF] Reed | PARTICIPATORY INDICATOR DEVELOPMENT: WHAT CAN ECOLOGISTS AND LOCAL COMMUNITIES LEARN FROM EACH OTHER[END_REF]. For this, robust indicators and methods that can be appropriated by stakeholders are needed. Such indicators should be relatively simple and thus rapid, robust and sensitive to the driver of interest, here climate [START_REF] Feld | The indication of biodiversity and ecosystem services: a synthesis across ecosystems and spatial scales[END_REF]. Sward height is a commonly used indicator of biomass production in ecological and agronomic studies [START_REF] Marriott | Long-term impacts of extensification of grassland management on biodiversity and productivity in upland areas[END_REF][START_REF] Oliveras | Grass allometry and estimation of above-ground biomass in tropical alpine tussock grasslands[END_REF][START_REF] Stewart | An evaluation of three quick methods commonly used to assess sward height in ecology[END_REF]. Accordingly, community mean plant vegetative height is a strong predictor of upland grassland production and its response to management [START_REF] Lavorel | Using plant functional traits to understand the landscape-scale distribution of multiple ecosystem services[END_REF][START_REF] Pakeman | Multivariate identification of plant functional response and effect traits in an agricultural landscape[END_REF]. Sward height is also regaining current interest given new imaging technologies that can capture fine-scale variations in grassland physical properties that can be linked to production and other ecosystem services [START_REF] Forsmoo | Drone-based structurefrom-motion photogrammetry captures grassland sward height variability[END_REF].

In this study we tested sward height as a method for monitoring mountain grassland production, with the objective that it should be useable for long-term monitoring by non-specialists (e.g. staff from protected areas, agricultural extension officers…). This paper first presents the sward heightbased field measurement method. We then test the robustness and sensitivity of this indicator by analysing a twelve-year time series representing a credible range of climate and productivity variation in the Central French Alps, for data collected across a set of 67 grassland plots from the main vegetation types across the Oisans / Briançonnais region [START_REF] Jouglet | Les végétations des alpages des Alpes françaises du Sud: guide technique pour la reconnaissance et la gestion des milieux pâturés d'altitude[END_REF]. Finally, we highlight how a potential tipping point in fodder production can be detected for recent summer temperature extremes. We conclude that desired criteria for an indicator of forage production response to climate variation are met and discuss possible extensions and improvements.

Materials and Methods

Field site

Plant biomass data was collected from a set of 67 grassland plots managed by farmers at the Lautaret study site (45°03'N, 6°24'E). The site is located in the Central French Alps on the southfacing slopes of Villar d'Arène and covers 13 km 2 with elevation ranging from1552 to 2442 m a.s.l. (a detailed site description can be found in [START_REF] Quétier | Plant traits in a state and transition framework as markers of ecosystem response to land-use change[END_REF]). It is managed by low-to medium intensity livestock rearing, involving a variable combination of organic fertilization at low doses (eight tons of manure per hectare every 2-3 years), mowing and grazing at low intensity (<2 days of livestock units per hectare per year). In total, six land-use types form a vegetation mosaic: three on previously cultivated terraces (currently fertilized and mown (TFM), mown but not fertilized (TM), or unmown and grazed in spring and autumn (TU)), two on permanent grasslands with no history of cultivation and a multi-century history of mowing (currently mown (UM), unmown and summergrazed (UU), and one on never mown summer grasslands (>2000 m) (SP). Average production at peak biomass (July) ranges from 2.5 to 5.5 tDM.ha -1 .yr -1 [START_REF] Jouglet | Les végétations des alpages des Alpes françaises du Sud: guide technique pour la reconnaissance et la gestion des milieux pâturés d'altitude[END_REF][START_REF] Lavorel | Using plant functional traits to understand the landscape-scale distribution of multiple ecosystem services[END_REF].

The climate is subalpine with a strong continental influence due to a rain shadow with respect to dominant westerly winds. Mean annual precipitation is 956 mm at the nearest Météo France climate station (Besse-en-Oisans, station number 38040001, Lat 6°10'18"E Long 45°04'12"N, Altitude 1525m) and the altitude-adjusted mean monthly temperatures of -4.6°C in January to 11°C in July (at 2050 m a.s.l.). Snow disappears from the study area between late March (at 1700 m) and mid-to late-May (at 2100 m). Regional mean temperatures have increased by 0.3°C per decade since the beginning of the 20 th century with a greater increase (0.4-0.5°C) in summer than in winter (0.1°C) or autumn (0.2°C) (GREC-Sud, 2018). Total annual precipitation is stable to slightly decreasing for 1959-2015, with a weak increase in spring precipitation. Overall the precipitation signal remains uncertain, as are precipitation projections. Interannual climate variability is marked and has increased since 2000 (Beniston, 2015;Durand et al., 2009b), including recurrence of heat waves [START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF].

Measuring standing biomass using sward height as an estimator

Reference biomass measurements

In the years 2007-2010 biomass estimates were carried out in a subset of twenty four representative 100 m² grassland plots using calibrated visual estimations of standing biomass [START_REF] Redjadj | Estimating herbaceous plant biomass in mountain grasslands: A comparative study using three different methods[END_REF]. In each plot, two operators made visual estimates of total (green biomass and litter) and green dry biomass in twelve independent randomly located 50x50 cm quadrats (i.e. 24 quadrats per plot across the two operators) each. These 24 quadrats were harvested manually to ground level, pooled and sorted to green and dead biomass, dried and weighed. For each operator during each year a calibration regression equation was then established between the operator's mean visual estimate of biomass per ha and the sampled standing and green biomass per ha. We note that this calibration was across plots and not at quadrat-level because individual quadrat visual estimates are not considered as meaningful on their own as opposed to plot averages [START_REF] Tothill | BOTANAL-a comprehensive sampling and computing procedure for estimating pasture yield and composition. 1. Field sampling[END_REF].The calibration equation was then used to correct the visual estimates in those quadrats in which actual biomass was not harvested (Suppl. Mat. A). Estimated biomass values for the 12 quadrats in each plot were averaged to give a mean total and green biomass per plot in t/ha per plot.

Relationship between reference biomass measurements and measures of sward height

Concurrently with these calibrated visual estimates of total and green biomass, sward height measurements were also taken within at least each of the 24 plots (total numbers 2007: 27, 2008: 41, 2009: 25, 2010: 26). Such measures of table or canopy closure height have been previously used in grasslands to successfully provide estimates of standing biomass [START_REF] Redjadj | Estimating herbaceous plant biomass in mountain grasslands: A comparative study using three different methods[END_REF][START_REF] Stewart | An evaluation of three quick methods commonly used to assess sward height in ecology[END_REF] and have been shown to reliably reflect spatial variations in biomass production across Swiss mountain grasslands [START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF]. Sward height measurements involved measuring 60-80 replicates of the height at which vegetation in the plot began to form a closed canopy (otherwise called table height). Measurements were averaged to give a mean sward height per plot. Measurements for the total of 119 plots across the four years (2007)(2008)(2009)(2010) showed a highly significant relationship between vegetation sward height and biomass (R² = 0.689 for total biomass, and 0.701 for green biomass), and no significant difference in this relationship between years (Suppl. Mat. A). Considering the six land-use types individually, while there was overall a highly significant relationship between sward height and biomass (R² = 0.697), and there was no significant difference in this relationship between land-use types (Suppl. Mat. A).

Validated method for using sward height to estimate total and green biomass Given the lack of significant differences across land-use types in the relationships between green biomass and sward height, data from all types and over the four years of measurements was pooled and a simple regression equation established explaining a high proportion of the variance between these two parameters (R 2 = 67%, p<0.001) (Figure 1). Measurements of biomass from 2010 onwards were carried out using only sward height measurements converted to green biomass using the following equation: Green Biomass (t DM/ha) = (0.0938 * sward height (cm)) + 1.2616 Figure 1. Reference calibration between peak green biomass and vegetation height (2007-2010)

Biomass data set

The data set consisted in such estimated biomass records from 2007 to 2018, representing 12 timecorrelated replicates. Mean green biomass was determined in 2-24 plots from each of six land-use types (TFM, TM, TU, UM, UU, SP) per year (due to plots being missed due to grazing/mowing, vole damage etc.), giving a total of 649 samples over the time course of the study, with 35-67 plots measured each year (Details in Suppl. Table B1). The plots in which this time-series biomass data was taken were the same as those in which the calibration equation was developed.

Climate data set

Climate data (daily minimum and maximum temperatures; precipitation) was obtained for 2007 to 2018 from the Météo France weather station at Besse-en-Oisans (Station number 38040001, Lat 6°10'18"E Long 45°04'12"N, Altitude 1525m) (Figure 2). This was the nearest complete station to the Lautaret field site and was assumed to provide an acceptable proxy for yearly climate variation. A temperature correction of +0.6°C per 100 m altitudinal gain was applied [START_REF] Theau | Herb'âge, une méthode de calcul des sommes de températures pour la gestion des prairies[END_REF] to inter-calibrate the station of Besse-en-Oisans (Altitude 1525 m) and a second station at Villar d'Arène Pied du Col (altitude 1665 m), and validated using a regression analysis over the 2007 growing season (1 April -31 July; R² = 0.662, n = 122). A further validation of Besse-en-Oisans as a reference station was made by a regression analysis between the Villar d'Arène meteorological station data for 2008 and records by two Hobo Onset® stations at the Les Cours locality (altitude 1810 m) and at the Lautaret alpine botanical garden (altitude 2070 m). Mean temperatures over the study period were ~2°C higher than the long term mean, with mean minimum monthly temperatures tracking most closely the long term average (Figure 2). There was a 2-3°C increase in spring and summer as compared to long term mean, while autumn and winter temperatures were less different. Mean precipitation was overall similar to the long-term average. A Principal Component Analysis of interannual variation across the study period is presented in Suppl. Mat. Figure C1. period (2007-18) as compared to the long term mean monthly temperatures . Vertical error bars represent the standard error of the mean monthly temperatures over the study period . B. Mean, maximum and minimum monthly precipitation over the study period at Besse-en-Oisans as compared to long term mean monthly precipitation . Vertical error bars represent the standard error of the mean monthly precipitation over the study period .

A.

B.

Synthetic climate variables were calculated in order to capture relevant agrometeorological parameters [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Nettier | Adaptation au changement climatique sur les alpages. Modéliser le système alpageexploitations pour renouveler les cadres d'analyse de la gestion des alpages par les sytèmes pastoraux[END_REF]. The cumulated growing degree days (GDD) over a time window is an indicator of cumulated energy for growth, and known to affect phenology of common grassland species (Ansquer et al., 2009). The cumulated number of frost events over the time window (FST; TN < 0°C) represents risks of damage to growing tissues, especially at the beginning of the season. Mean monthly soil moisture (MSMT) was estimated using the Thornthwaite monthly water balance model [START_REF] Garnier | A standardized methodology to assess the effects of land use change on plant traits, communities and ecosystem functioning in grasslands[END_REF][START_REF] Thornthwaite | The water balance[END_REF]. A remote sensing estimate of the date of snowmelt each year at site level was added as a proxy for of the date of onset of the growing season. Chosen climatic explanatory variables thus encompassed a range of known potential climatic drivers of differences in biomass production in alpine herbaceous vegetation, including variables capturing temperature, soil moisture, frost and growing period length (snow free period), and this over time windows relevant to determining plant growth and biomass development in a sub-alpine environment [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF]. Each variable was calculated for differing time windows of each year so as to isolate parts of the growth period considered a priori most important in determining plant growth and biomass development. Suppl. Mat. C presents full list of calculated parameters (Table C1), with an explanation of the calculation methods and definition of the various time windows. Initial analyses using climate variables calculated over the whole year, or for the full seasons of Winter (Dec, Jan, Feb), Summer (June, July, August) or Autumn (Sep, Oct, Nov) showed no significant correlations with the biomass variables, and thus were not included in the final analyses. Likewise, an exploration of the relevance of decadal rather than monthly resolution for explanatory climate variables showed no improvement on statistical models, hence with operationality for users we kept the monthly resolution.

Analysis of correlations between differences in biomass production and the climatic variables:

For the correlations with climate variables, as all of the same plots were not necessarily measured every year, and most of all because there was only one climate measure common to all plots each year, it was not possible to use individual plots as independent data points. For each land-use type the mean biomass of all plots measured in a given year was calculated, and then these values averaged across the 11 years from 2007 to 2017 to give a mean biomass per land-use type over the time course of the survey. 2018 was only incorporated in final model validation (see below). Mean biomass per land-use type for the whole time period was then subtracted from the mean biomass per land-use type of each year, to give the biomass production anomaly from the overall mean for each year. Thus the biomass data used in the climate correlations consisted of a value of the mean biomass anomalies from the 11-year average mean biomass per trajectory per year (thus 11 annual data points per land-use type). Two other biomass variables were also calculated: mean biomass anomaly per year for grasslands other than summer pastures (TFM, TM, TU, UU and UM combined; all < 2000 m), and the mean biomass anomaly per year for all types combined.

The data set thus comprised of series of 8 response variables (the biomass anomalies per year for the six land-use types, plus the two combined biomass anomalies), and a series of 40 explanatory variables, over the 11 years of measurement. To determine what part of the variation in anomalies in biomass production between years could be attributed to climatic effects, and which climate variables were explanatory, we used linear regression with correlated errors using REML estimation (Genstat 16 th edition). This method can be used to fit regression models to data such as repeated measurements. Where required data was transformed to meet assumptions of normality. For each biomass response variable we initially carried out a simple regression with each explanatory climate variable in order to select relevant variables (Suppl. Table 1). Those climate variables which had a significant correlation with biomass anomalies were retained and ranked in order of their degree of variance explained. Subsequently, multiple regression models for each of the biomass anomaly variables were developed by fitting those climate variables significant on their own in a stepwise regression manner, adding variables in the order of their degree of variance explained singly. Added variables were retained, or discarded on the basis of their significance in the presence of other, more significant variables, until a final most parsimonious model on the basis of compared Aikake criterion values was retained which included only significant climate variables and which maximized the total variance explained.

In order to synthesise effects and their variations across land-use types a co-inertia analysis of the biomass variation over time amongst the land-use types (mean biomass anomaly each year by type) with the climate variation over time (a selection of climate variables identified as being significant in the regression analyses) was carried out using the ade4 R-package software [START_REF] Dray | The ade4 package: implementing the duality diagram for ecologists[END_REF].

Results

Interannual variability in peak biomass

Mean peak biomass production varied from an average of 1.5 -3.25 t/ha (kg/m²) for summer pastures (SP) and grazed terraces (TU), to 4.1 -6.6 t/ha (kg/m²) for fertilized hay meadows on terraces (TFM) (Figure 3). Coefficients of variation were greatest for TFM (44%) and UM (31%) and smallest for TM (21%) (Suppl. Mat. B, Figure B2). In the following we first report the explanatory power of statistical models of inter-annual variation in peak biomass and the relative contributions of different significant variables for successive analyses of the complete data set (all plots), all land-use types except summer pastures (due to their qualitatively different response) and individual land-use types. A first highly explanatory model was developed for years until 2014, and is henceforth considered as the baseline period. Inclusion of summer months (June-July) for subsequent years, including the three warmest years on record (2015,2017,2018) resulted in no significant models being found. The final models retained for the baseline period (Table 1) revealed the strong effects of air temperatures during the period of maximum plant growth (June-July) in determining peak biomass. Mean maximum temperatures during the period of maximum plant growth (June-JulyTX) explained 61.5% of inter-annual variance across all plots, as compared to 49.1% for growing degree days for the same period (June-JulyGDD) (Suppl. Table 1). When excluding summer pastures, mean maximum temperatures during July (JulyTX) explained 60.6% of inter-annual variance , as compared to 49.1% for growing degree days for the same period (July GDD). The final best fit multi-linear models explained high amounts of variation in the biomass production, with mean minimum temperature across the growth period and July precipitation explaining 89% of inter-annual variance across all plots. In models for individual landuse types the percentage of variance explained ranged from 36% (UU) to 92% (TFM), and reached 86% for grasslands other than summer pastures.

Response of biomass production anomalies to climate variables

Secondly, we detail the effects of significant climate variables. In the analysis considering land-use types other than summer pastures, there was a significant positive effect on peak biomass of increasing temperatures (TN or TX) over the growth period until peak biomass and especially in June-July, the most active plant growth period (Figure 4, Suppl. Fig. 1). Statistical models for individual land-use types retained either mean minimum (TFM, TM, grasslands other than summer pastures) or mean maximum (TM, TU, UM, UU) temperatures (Table 1). Growing degree days were never retained in final models given their less explanatory effects than temperature variables (Suppl. Table 1), although for TFM a model including JulyGDD was nearly as explanatory as the final model with mean minimum temperatures (Table 1). We also observed that increase in precipitation during this same growth period was a weak predictor of decreases in peak biomass in TFM and TU. This effect is most likely due to significant correlations (R² = 0.453, p = 0.047 for Spring TX vs Spring RR during the baseline period) between increasing precipitation and cooler low-radiation conditions (Suppl. Mat. C Figure C2), resulting in decreasing plant growth and biomass production beyond lower temperatures effects. Increasing frosts during May significantly reduced peak biomass in UM and especially SP located at higher altitude.

Co-inertia analysis

The co-inertia analysis of the biomass variation over time amongst land-use types (mean biomass anomaly each year by land-use type) with the variation over time of significant climate variables showed consistent patterns with the regressions. The RV-coefficient showed a highly significant correlation of 0.744 between the two analyzed tables (Monte-Carlo test based on 999 replicates, simulated p-value: 0.002). The first eigenvalue of the co-inertia analysis explained 85.8 % of the total explained variance between the two datasets, showing primary relevance of the links between biomass variation and the climate parameters along this first axis.

This analysis emphasized the positive associations revealed by the univariate analyses between peak biomass anomalies (especially in types TFM, TM, TU and UM) with growth period and June-July mean maximum temperatures, and the negative association of anomalies with increasing accumulated precipitation over the growth period or June-July (Figure 5). Also, the analysis singled out the response of summer pastures whose negative biomass anomalies were primarily determined by May frosts. Note also the lack of responsiveness of unmown and summer-grazed (UU) grasslands to analysed variables, consistent with univariate analyses (Suppl. Table 1). C1.

Anomalies in the hottest summers

We failed to find any significant model when inserting one or all of the following years: 2015, 2017 and 2018, which were the hottest summers within our time series (mean June-July TX of 23.8, 22.8 and 22.5°C respectively). Their deviation from the baseline relationships with temperature could therefore not be captured by any of the available variables. Adding them as supplementary points in the baseline regression models (Figure 4, Suppl. Fig. 1 for complete results) shows that for the three years where the mean June-July maximum temperatures exceeded 20°C, and which were clearly anomalous from the baseline regime in terms of weather (Suppl. Figure C1), the positive effect of increasing temperature was lost.

Discussion

Climatic controls of biomass production anomalies

Our analyses of responses of peak biomass to climate variation across 12 years underline the dominant effects of temperature during the peak growth period (June and July) in determining the production of biomass in mountain grasslands. Given the coarseness of the analysis which used biomass estimated via the sward height proxy, climate data from a distant (yet representative) site and land-use type mean data for biomass production, the significance and degree of variance explained of these models is remarkable and supports the robustness and sensitivity of sward height as an indicator for climate response of peak fodder biomass.

Sums of temperature (or growing degree days -GDD) are frequently used for capturing, possibly in combination with radiation or its indicators (modifications by slope and aspect - [START_REF] Jäger | Mapping of supply and demand of biomass from alpine grassland[END_REF]), available energy that all physiological processes, and thus growth depend on. Plant phenology is also strongly linked with thresholds in GDD, thereby determining critical periods of vegetative growth and reproduction [START_REF] Dumont | A metaanalysis of climate change effects on forage quality in grasslands: specificities of mountain and Mediterranean areas[END_REF]Duru et al., 2010b;[START_REF] Henebry | Phenologies of North American Grasslands and Grasses[END_REF]. [START_REF] Jonas | How alpine plant growth is linked to snow cover and climate variability[END_REF] found a positive effect of GDD on biomass production through the snow-free season (considered as the entire growth period) and a more moderate effect of temperature during the month following snow melt across 17 alpine sites in Switzerland, but their analysis did not differentiate months within the growth season. In general, where there were significant relationships between TN or TX and biomass difference there was also a significant relationship with GDD, but it was weaker than respective TN / TX relationships (Suppl. Table 1). Overall, average growth period temperatures were more relevant to TM, TU, UM and across land-use types (Table 1, Figure 5). But given the strong correlation between temperature variables (Suppl. Mat. C, Figure C1), the overall thermal energy response may be captured by either. Although GDD may be more integrative, given response thresholds identified with TX (see below) we nevertheless suggest using them in combinations until specific effects and mechanisms are confirmed.

In our analyses snowmelt date (ONSET) was never retained as a significant variable, contrary to [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF] at regional level for analyses of growth season NPP. Snowmelt date did actually vary by one month across years (mean snow melt date at site level from the 24 th April in 2011 to the 23 rd May in 2013) and was only weakly related to peak biomass for higher pastures (SP and UU) (Suppl. Table 1). Yet, snowmelt date was well correlated with mean minimum temperature in May (R² = 0.657) and associated number of frosts in May (R² = 0.697), which, consistent with expectations (Liu et al., 2018b), decreased peak biomass in less productive vegetation types (TU, UM, SP; Figure 5). We suspect that the decadal extent of our data set was not enough to clearly detect a direct significant response to snowmelt date that may be captured by longer (e.g. 30 year) time series [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF]. Additionally, soil rather than air temperature may be a more relevant variable to alpine plant responses to climate [START_REF] Choler | Winter soil temperature dependence of alpine plant distribution: Implications for anticipating vegetation changes under a warming climate[END_REF][START_REF] Guo | Early-spring soil warming partially offsets the enhancement of alpine grassland aboveground productivity induced by warmer growing seasons on the Qinghai-Tibetan Plateau[END_REF], but requires onsite equipment that may not be widely available. Our results nevertheless suggest that simple meteorological variables such as mean temperatures may still be sufficiently informative to account for peak biomass responses.

While drought is considered as a main limiting factor for temperate and mountain grassland production (Calanca et al., 2016;[START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF]Duru et al., 2010a;[START_REF] Trnka | A simple statistical model for predicting herbage production from permanent grassland[END_REF], our analyses did not identify any negative effects of reduced precipitation or soil moisture during the growth period to peak biomass. Soil water recharge from snowmelt may support the first growth cycle regardless of precipitation. We however note that our time series did not include any radically dry springs. The three driest years 2009, 2010, 2012 in our data set did not have reductions in production as compared to prediction by June-July temperatures. The future recurrence of very low snow years (such as 2017) [START_REF] Verfaillie | Multi-component ensembles of future meteorological and natural snow conditions for 1500 m altitude in the Chartreuse mountain range, Northern French Alps[END_REF] may however challenge this buffering capacity. Additionally, drought is likely to affect post-harvest or grazing regrowth [START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF][START_REF] Sautier | Assessing climatic exposure of grasslandbased livestock systems with seasonal-scale indicators[END_REF], especially in the context of extended growing seasons by increasing later summer temperatures [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF], but such effects could not be captured with our assessment of climate effects on peak biomass.

Lastly, different land-use types differed in their specific sensitivities to climate parameters. Temperature effects (effect sizes in Table 1) were strongest in TFM, TM and SP, where dominant plants with more exploitative strategies (e.g. lower Leaf Dry Matter Content -LDMC; [START_REF] Quétier | Plant traits in a state and transition framework as markers of ecosystem response to land-use change[END_REF]) may be expected to have more sensitive growth to environmental variation such as temperature and its effects on soil nutrient and moisture availability [START_REF] Jung | Intraspecific trait variability mediates the response of subalpine grassland communities to extreme drought events[END_REF][START_REF] Liancourt | Leaf-trait plasticity and species vulnerability to climate change in a Mongolian steppe[END_REF][START_REF] Stampfli | High land-use intensity exacerbates shifts in grassland vegetation composition after severe experimental drought[END_REF]. Conversely, the more conservative dominant plants (higher LDMC) with lower phenotypic plasticity found in TU and UU [START_REF] Grassein | Plant resource use strategies: the importance of phenotypic plasticity in response to a productivity gradient for two subalpine species[END_REF] may explain the lower temperature responsiveness of peak biomass in these grasslands [START_REF] Karlowsky | Land use alters drought responses of carbon allocation and plant-microbial interactions in mountain grassland[END_REF]. Also differences across land-use types in climate specific parameter sensitivities mean that different grasslands have asynchronous good and bad years. These differential sensitivities are important for farmers adaptation to climate variability because they imply that different land-use types may be used as resources and adjustments within and across years (Andrieu et al., 2007;[START_REF] Nettier | Resilience as a framework for analysing the adaptation of mountain summer pasture systems to climate change[END_REF]. Specifically, hay from a warm year in mown terraces can buffer e.g. a subsequent poor spring; extreme early summer heat may affect lower grasslands including hay meadows but not summer pastures (if not associated with drought - [START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF]); and insensitive grasslands (UU) offer a buffer when a cold season reduces production in other types.

Anomalies in the three hottest summers on record

Peak biomass production deviated strongly from the temperature-driven relationship for the three hottest summers in our data set (Figure 4). Other variables, especially early season conditions (May temperatures, frosts) showed responses aligned with other years (Suppl. Fig. 1), but no alternative model incorporating these three years could be found. This suggests the possibility of a regime shift, where baseline driver responses are departed from [START_REF] Folke | Regime Shifts, Resilience, and Biodiversity in Ecosystem Management[END_REF]. Such departures from the baseline regime would obviously need confirmation across longer time series and more sites. Nevertheless, they may be related to known grassland growth mechanisms.

Firstly there are known heat response thresholds for cool temperate grasslands from 18°C for grasses to 23°C for legumes (Duru et al., 2010a), and an overall threshold of 20°C has been considered for growth modelling (Calanca et al., 2016). [START_REF] Cremonese | Heat wave hinders green wave: The impact of climate extreme on the phenology of a mountain grassland[END_REF] indeed observed severe reduction in photosynthetic activity and growth for Nardus stricta grassland during the summer 2015 heat wave in the Italian Alps. Additionally physiological damage from oxidative stress at peak temperatures in the sun, which can excess 40°C in alpine grasslands [START_REF] Körner | Alpine plant life. Plant ecology of high mountain ecosystems[END_REF], has been shown to occur for cool temperate tussock grasses (e.g. Lolium perenne, [START_REF] Soliman | Oxidative stress and physiological damage under prolonged heat stress in C3 grass Lolium perenne[END_REF]). Secondly, heat accelerates phenological development [START_REF] Dumont | A metaanalysis of climate change effects on forage quality in grasslands: specificities of mountain and Mediterranean areas[END_REF], and may thereby shorten the effective growth period from snow melt, with possible repercussion on maximum vegetative growth. While physiological limitations are likely, we consider this as the most parsimonious explanation for the growth shortage observed during the three hottest years. At our site, we previously observed flowering for the most common grasses to be temperature-driven (e.g. 638 °C.days for dominants Dactylis glomerata and Bromus erectus in terraced grasslands TFM, TM, TU; 494 °C.days for the dominant Patzkea paniculata in UM and UU grasslands; Lavorel et al. unpublished). In mountain grasslands, the active vegetative growth period is considered, consistent with these phenological thresholds, to be on average between 300 and 600 °C.days [START_REF] Nettier | Adaptation au changement climatique sur les alpages. Modéliser le système alpageexploitations pour renouveler les cadres d'analyse de la gestion des alpages par les sytèmes pastoraux[END_REF]. Its duration has been continuously decreasing over the years in our data set, with 2015 and 2017 being the two shortest (Suppl. Mat. C, Figure C3).

Further, [START_REF] De Boeck | Simulated heat waves affected alpine grassland only in combination with drought[END_REF] found no experimental effect of heat per se on growth of alpine vegetation turfs unless it was associated with drought, consistent with other temperate grasslands [START_REF] Hoover | Resistance and resilience of a grassland ecosystem to climate extremes[END_REF]). An analysis across the entire French Alps confirmed that heat waves were associated with reduced annual grassland production only for drought years, and specifically due to lack of summer regrowth [START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF]. At our site, in 2017 and 2018 growing seasons until July were relatively dry as compared to the average of the study period, while 2015 was rather wet (e.g. July RR, Suppl. Mat. C, Figure C1A). Nevertheless, the precipitation-related deviation in peak biomass did not align with other years even in fertilised hay meadows (TFM; most sensitive to drought - [START_REF] Karlowsky | Land use alters drought responses of carbon allocation and plant-microbial interactions in mountain grassland[END_REF]) (negative regression with July RR, Supplementary Figure 1). We conclude that the three recent extreme years likely each represent unique combinations of direct heat effects, shortened vegetative growth season and drought, but that their recurrence raises concern for a climate tipping point in mountain grassland fodder production. The baseline regime and the anomalies in hottest years need confirmation for longer time series and across a greater geographic extent.

Implications for fodder production, ecosystem services and farming in mountains

Our analyses showed significant effects of growing season temperatures on mountain fodder production, with differences in driving parameters and sensitivities across land-use types. We also revealed the potential for extreme summer temperature to reduce potential gains from increased average temperatures and decreased spring frosts. Recent observed trends and projections for these climate parameters across the Alps and other European mountains (Beniston et al., 2018;Durand et al., 2009a;[START_REF] Hock | IPCC Special Report on Ocean and Cryosphere in a Changing Climate[END_REF][START_REF] Zubler | Key climate indices in Switzerland; expected changes in a future climate[END_REF] thus suggest potential gains in fodder production from increased mean spring temperatures (from +2 to +6°C), which will trade-off with risks from spring frosts and more frequent summer temperature extremes as those observed in the most recent years [START_REF] Corona-Lozada | Drought offsets the positive effect of summer heat waves on the canopy greenness of mountain grasslands[END_REF][START_REF] Cremonese | Heat wave hinders green wave: The impact of climate extreme on the phenology of a mountain grassland[END_REF]. Risk from spring frost is of particular concern especially for sensitive high elevation summer pastures given declining length of the snow season, meaning an early start in vegetation growth and potential exposure to damaging frost [START_REF] Choler | Growth response of temperate mountain grasslands to inter-annual variations in snow cover duration[END_REF][START_REF] Klein | Unchanged risk of frost exposure for subalpine and alpine plants after snowmelt in Switzerland despite climate warming[END_REF]. Recurrence of extreme summer temperatures induces further uncertainties, with effects yet to be confirmed across the Alps, along with their magnitude.

Nevertheless, we emphasized the potential for resilience for farms from their multiple land-uses across the landscape, allowing them to balance shortages from some grasslands with stability or gains from others depending on years. Farmers and herders also use a diversity of management strategies within seasons, years and across years supporting their resilience, whether from adjustment in herd sizes, grazing practices (e.g. adjusting grazing dates and spatial distribution) or from stock management and potential purchase of hay [START_REF] Nettier | Resilience as a framework for analysing the adaptation of mountain summer pasture systems to climate change[END_REF]. Additionally, because the response of biomass production to climate change underpins those of multiple grassland ecosystem processes and services like soil carbon sequestration, regulation of water quantity and quality or habitat for invertebrates [START_REF] Lavorel | Climate change effects on grassland ecosystem services[END_REF], observed trends are likely to cascade to multiple impacts for farmers and other land managers in mountains. We conclude that long-term monitoring of mountain grassland biomass production with simple methods and indicators, and analyses of its response to climate parameters is critical for supporting farmers and land managers climate change adaptation [START_REF] Deléglise | A Method for Diagnosing Summer Mountain Pastures' Vulnerability to Climate Change, Developed in the French Alps[END_REF].

Acknowledgements

This research was carried out as part of LTSER Zone Atelier Alpes and the Sentinelles des Alpes -Alpages Sentinelles Programme, with joint funding from French Agence Nationale pour la Recherche (ANR) MtnPaths (ANR-16-CE93-0008-01). This publication contributes to ANR Investissements d'Avenir CDP Trajectories (ANR-15-IDEX-02).

List of Tables and Figures

Table 1. Climate variables retained within the multiple regression models explaining peak biomass variation for each of the land-use types. For each of them the % variation explained by the retained climate variables, the significance (P) of the full model, the significance of the effect of each retained climate variable and the direction and magnitude of the standardised effect for each of the retained climate variables are presented. TN: mean minimum temperature, TX: mean maximum temperature, RR: cumulated precipitation, FST: number of frosts, GrthPrd: growth period, SMST: soil moisture index.

Andrieu, N., Josien, E., Duru, M., 2007

Figure 2 .

 2 Figure 2. A. Mean, maximum and minimum monthly temperatures at Besse-en-Oisans over the studyperiod (2007-18) as compared to the long term mean monthly temperatures. Vertical error bars represent the standard error of the mean monthly temperatures over the study period(2007- 18). B. Mean, maximum and minimum monthly precipitation over the study period at Besse-en-Oisans as compared to long term mean monthly precipitation. Vertical error bars represent the standard error of the mean monthly precipitation over the study period.

Figure 3 -

 3 Figure 3 -Time series of mean peak biomass across land-use types. Error bars represent variation across plots for a given land-use type.

Figure 5 .

 5 Figure 5. Correlations between selected climate variables (vertical) and peak biomass anomalies across land-use types (horizontal) from co-inertia analysis. The size of individual squares represents the strength of the correlation. Abbreviations for climate variable are listed in TableC1.

Figure 1 .

 1 Figure1. Reference calibration between peak green biomass and vegetation height(2007)(2008)(2009)(2010)

Figure 3 .

 3 Figure 3. Time series of mean peak biomass across land-use types. Error bars represent variation across plots for a given land-use type.

Figure 4 .

 4 Figure 4. Relationships between peak biomass and significant explanatory variables for A. All landuse types, B. All grasslands except summer pastures, C. Summer pastures. Each graph presents data with significant regression for the baseline regime (2007-2014+2016) and additional points for heat extremes of 2015, 2017 and 2018. The full set of relationships for individual land-use types is presented in Supplementary Figure 1.

Figure 5 .

 5 Figure 5. Correlations between selected climate variables (vertical) and peak biomass anomalies across land-use types (horizontal) from co-inertia analysis.

Table 1 .

 1 Climate variables retained within the multiple regression models explaining peak biomass variation for each of the land-use types. For each of them the % variation explained by the retained climate variables, the significance (P) of the full model, the significance of the effect of each retained climate variable and the direction and magnitude of the standardised effect for each of the retained climate variables are presented. TN: mean minimum temperature, TX: mean maximum temperature, RR: cumulated precipitation, FST: number of frosts, GrthPrd: growth period, SMST: soil moisture index.

	Biomass	Retained	% variation	Full model	Effects p	Standard.
	anomaly	climate	explained by full	p		Effect
	response	variables	model
	variable
	TFM	JulyTN	91.6	<.001	0.005	0.4140
	(TN model)	JulyRR			0.043	-0.00743
	TFM	JulyGDD	90.1	<.001	0.009	0.01165
	(GDD model)	JulyRR			0.039	-0.00592
	TM	GrthPrdTN	45.8	0.027	0.027	0.444
	(TN model)
	TM	GrthPrdTX	37.7	0.046	0.046	0.2176
	(TX model)
	TU	GrthPrdTX	68.5	0.013	0.044	0.270
		JuneJulyRR			0.035	-0.00537
	UM	JuneJulyTX	79.7	0.004	0.002	0.2974
		MayFST			0.016	-0.0764
	UU	GrthPrdTX	35.6	0.053	0.053	0.474
	SP	MayFST	42.1	0.035	0.035	-0.0870
	T1-5	JulyTX	60.6	0.008	0.008	0.1457
	(Simple model)
	T1-5	GrthPrdTN	85.7	0.001	0.013	0.412
	(Complex	JulySMST			0.003	-0.00881
	model)
	All	JuneJulyTX	61.5	0.008	0.008	0.2400
	(Simple model)
	All	GrthPrdTN	89.2	<.001	0.052	0.2251
	(Complex	JulyRR			0.001	-0.005570
	model)
	Figure 4. Relationships between peak biomass and significant explanatory variables for A. All land-
	use types, B. All grasslands except summer pastures, C. Summer pastures. Each graph presents data
	with significant regression for the baseline regime (2007-2014+2016) and additional points for heat
	extremes of 2015, 2017 and 2018. The full set of relationships for individual land-use types is

presented in Supplementary Figure

1

. Abbreviations for climate variable are listed in Table

C1

.

  . Relationships between diversity of grassland vegetation, field characteristics and land use management practices assessed at the farm level. Agriculture, Ecosystems and Environment 120, 359-369. Ansquer, P., Al Haj Kahled, R., Cruz, P., Theau, J.P., Thérond, O.,Duru, M., 2009. Characterizing and predicting plant phenology in species rich grasslands. Grass and Forage Science 64, 57-70. Beniston, M., 2015. Ratios of record high to record low temperatures in Europe exhibit sharp increases since 2000 despite a slowdown in the rise of mean temperatures. Climatic Change 129, 225-237. Beniston, M., Farinotti, D., Stoffel, M., Andreassen, L.M., Coppola, E., Eckert, N., Fantini, A., Giacona, F., Hauck, C., Huss, M., 2018. The European mountain cryosphere: a review of its current state, trends, and future challenges. Cryosphere 12, 759-794. Calanca, P.,Deléglise, C., Martin, R., Carrère, P., Mosimann, E., 2016. Testing the ability of a simple grassland model to simulate the seasonal effects of drought on herbage growth. Field Crops Research 187, 12-23. Carlson, B.Z., Choler, P., Renaud, J., Dedieu, J.-P., Thuiller, W., 2015. Modelling snow cover duration improves predictions of functional and taxonomic diversity for alpine plant communities. Annals ofBotany 116, 1023Botany 116, -1034. .