Iñigo Gómara
email: i.gomara@upm.es

Gianni Bellocchi

Raphaël Martin

Belén Rodríguez-Fonseca

Margarita Ruiz-Ramos

Influence of climate variability on the potential forage production of a mown permanent grassland in the French Massif Central

Climate Services (CS) provide support to decision makers across socio-economic sectors. In the agricultural sector, one of the most important CS applications is to provide timely and accurate yield forecasts based on climate prediction.

In this study, the Pasture Simulation model (PaSim) was used to simulate, for the period 1959-2015, the forage production of a mown grassland system (Laqueuille, Massif Central of France) under different management conditions, with meteorological inputs extracted from the SAFRAN atmospheric database. The aim was to generate purely climate-dependent timeseries of optimal forage production, a variable that was maximized by brighter and warmer weather conditions at the grassland.

A long-term increase was observed in simulated forage yield, with the 1995-2015 average being 29% higher than the 1959-1979 average. Such increase seems consistent with observed rising trends in temperature and CO 2 , and multi-decadal changes in incident solar radiation. At interannual timescales, sea surface temperature anomalies of the Mediterranean (MED), Tropical North Atlantic (TNA), equatorial Pacific (El Niño Southern Oscillation) and the North Atlantic Oscillation (NAO) index were found robustly correlated with annual forage yield values. Relying only on climatic predictors, we developed a stepwise statistical multi-regression model with leave-one-out cross-validation. Under specific management conditions (e.g., three annual cuts) and from one to five months in advance, the generated model successfully provided a p-value < 0.01 in correlation (t-test), a root mean square error percentage (%RMSE) of 14.6% and a 71.43% hit rate predicting above/below average years in terms of forage yield collection. This is the first modeling study on the possible role of large-scale oceanic-atmospheric teleconnections in driving forage production in Europe. As such, it provides a useful springboard to implement a grassland seasonal forecasting system in this continent.

Introduction

Nowadays, a great demand exists for translating the large amount of climate data and information into customized tools, products and services (European Commission, 2014). These tools, so-called Climate Services (CS), provide support to decision makers across socio-economic sectors (e.g., water management, agriculture, energy production). In the agricultural sector, one of the most important CS applications is to provide timely and accurate yield forecasts based on climate prediction, which can be used to adopt best management practices and readily incorporate them into on-farm management plans in a cost-effective manner [START_REF] Lemos | Narrowing the climate information usability gap[END_REF][START_REF] Mase | Unrealized potential: a review of perceptions and use of weather and climate information in agricultural decision making[END_REF]. User-driven climate-information services may also serve to prevent emergencies in target food-insecure regions (e.g. [START_REF] Haile | Weather patterns, food security and humanitarian response in sub-Saharan Africa[END_REF] or orient market forces and regulations related to food production (e.g. [START_REF] Whitfield | Frontiers in climate smart food systems: outlining the research space[END_REF], and support environmental protection while adapting to climate change (e.g. [START_REF] Hackenbruch | Tailoring climate parameters to information needs for local adaptation to climate change[END_REF].

Accurate weather predictions beyond a two-week horizon are not possible because of the intrinsic chaotic behavior of the Earth's atmospheric system [START_REF] Lorenz | Deterministic non periodic flow[END_REF][START_REF] Lorenz | Atmospheric predictability experiments with a large numerical model[END_REF][START_REF] Yoden | Atmospheric predictability[END_REF]. However, seasonal predictions can be retrieved using as input oscillatory and slowly https://doi.org/10.1016/j.agrformet. 2019.107768 varying components in the coupled ocean-atmosphere-land system that dominate the weather response [START_REF] Luo | Assessing the idealized predictability of precipitation and temperature in the NCEP climate forecast system[END_REF]. Sea surface temperatures (SSTs) are likely the most important source of climate predictability at multiple timescales (e.g., seasonal, interannual, multidecadal) because ocean water has a much higher heat capacity than air [START_REF] Bjerknes | Atmospheric teleconnections from the equatorial pacific[END_REF]. SST anomalies over large marine areas can affect the atmospheric circulation, so that some regions may become warmer and wetter than normal, while others become cooler and drier during specific periods [START_REF] Hurrell | Decadal climate prediction: challenges and opportunities[END_REF]. SST anomalies thus provide valuable information on the evolution of air-sea interactions and associated impacts in climate, in both nearby and distant locations [START_REF] Hand | The role of local sea surface temperature pattern changes in shaping climate change in the North Atlantic sector[END_REF], as these anomalies can alter and propagate into oceanic-atmospheric pathways (e.g. [START_REF] Wallace | Teleconnections in the geopotential height field during the northern hemisphere winter[END_REF][START_REF] Årthun | Skillful prediction of northern climate provided by the ocean[END_REF][START_REF] Gleixner | Seasonal predictability of Kiremt rainfall in coupled general circulation models[END_REF][START_REF] Mamalakis | A new interhemispheric teleconnection increases predictability of winter precipitation in southwestern US[END_REF]. This propagation takes place through atmospheric large-scale recurrent and persistent anomaly patterns, which are called teleconnection patterns (hereafter referred as "teleconnections" for simplicity; [START_REF] Barnston | Classification, seasonality and persistence of lowfrequency atmospheric circulation patterns[END_REF].

In the Euro-Atlantic sector, the North Atlantic Oscillation (NAO) is the most prominent teleconnection affecting climate variability on a broad frequency band (daily to multi-decadal). The NAO characterizes the mean sea level pressure (mslp) gradient between sub-polar and subtropical latitudes and is driven by multiple variability sources (like stochastic internal, oceanic, solar, stratospheric; [START_REF] Hurrell | The North Atlantic Oscillation: climate significance and environmental impact[END_REF][START_REF] Pinto | Past and recent changes in the North Atlantic Oscillation[END_REF]. Under its positive phase (NAO+), synoptic weather systems are deflected northward, causing cloudier, wetter conditions over northern Europe and sunnier, drier conditions over southern Europe [START_REF] Gómara | On the relation between explosive cyclones affecting Europe and the North Atlantic Oscillation[END_REF]. The reverse holds for NAO-.

The warming or cooling of Pacific Equatorial SSTs, known as El Niño Southern Oscillation (ENSO), is the major source of climate variability affecting different regions of the world (e.g. [START_REF] Trenberth | The definition of El Niño[END_REF]. At interannual timescales, a positive ENSO phase tends to be associated with NAO-in late winter [START_REF] Brönnimann | Impact of El Niño-Southern Oscillation on European climate[END_REF], and with increased precipitation over southern Europe during summer and fall [START_REF] Shaman | The seasonal effects of ENSO on european precipitation: observational analysis[END_REF], although this link appears to be non-stationary in time [START_REF] López-Parages | Multidecadal modulation of EL Niño influence on the Euro-Mediterranean rainfall[END_REF]. Additional SST patterns with influence over Europe are the Tropical North Atlantic (TNA) and the Mediterranean (MED). A warmer TNA is known to increase precipitation over southwestern Europe in winter [START_REF] Okumura | Tropical Atlantic air-sea interaction and its influence on the NAO[END_REF][START_REF] Czaja | Observed impact of atlantic SST anomalies on the North Atlantic Oscillation[END_REF][START_REF] Rodríguez-Fonseca | Evaluation of the North Atlantic SST forcing on the European and North African winter climate[END_REF][START_REF] Losada | Impacts of SST anomalies on the North Atlantic atmospheric circulation: a case study for the northern winter 1995/1996[END_REF]. A warmer MED is linked to higher temperatures over Europe in summer and autumn [START_REF] Feudale | Role of Mediterranean SST in enhancing the European heat wave of summer 2003[END_REF][START_REF] García-Serrano | Large-scale atmospheric response to eastern Mediterranean summer-autumn SST anomalies and the associated regional impact[END_REF]. At lower frequencies, the Atlantic Multidecadal Oscillation (AMO; a multidecadal warming/ cooling of North Atlantic SSTs) is known to influence the NAO [START_REF] Gastineau | Influence of the North Atlantic SST variability on the atmospheric circulation during the twentieth century[END_REF], the north Atlantic storm track [START_REF] Woollings | Response of the North Atlantic storm track to climate change shaped by ocean-atmosphere coupling[END_REF][START_REF] Gómara | Abrupt transitions in the NAO control of explosive North Atlantic cyclone development[END_REF][START_REF] Diodato | Enhanced propagation of rainfall kinetic energy in the UK[END_REF][START_REF] Vaideanu | Atlantic multidecadal oscillation footprint on global high cloud cover[END_REF] and Euro-Mediterranean precipitation [START_REF] López-Parages | Multidecadal modulation of EL Niño influence on the Euro-Mediterranean rainfall[END_REF][START_REF] Diodato | Using historical precipitation patterns to forecast daily extremes of rainfall for the coming decades in Naples (Italy)[END_REF]. Thus, teleconnections can provide prediction of key variables for agro-ecosystems like maximum/minimum daily temperatures, incident solar radiation, precipitation or wind speed, which can be assimilated into monthly-to-decadal forecasting systems (e.g. [START_REF] Dannenberg | Atmospheric teleconnection influence on North American land surface phenology[END_REF].

Decades of research have disclosed strong connections between ENSO and major crop yields worldwide. [START_REF] Cane | Forecasting Zimbabwean maize yield using eastern equatorial pacific sea surface temperature[END_REF] found a strong anticorrelation between ENSO index and maize yields in Zimbabwe (even stronger than for rainfall). [START_REF] Phillips | Recent advances in understanding the role of the El Niño-Southern Oscillation (ENSO) in climate variability present opportunities for improving efficiency in agricultural production[END_REF] found a positive effect of El Niño episodes in the U.S. corn belt, a result subsequently corroborated by [START_REF] Kellner | Climate variability and the U.S. corn belt: ENSO and AO episode-dependent hydroclimatic feedbacks to corn production at regional and local scales[END_REF]. [START_REF] Meinke | Using seasonal climate forecasts to manage dryland crops in northern Australia -experiences from the 1997/98 seasons[END_REF] showed that positive ENSO events depressed rainfall in Eastern Australia, leading to lower-than-average wheat yields. Other SST indices were investigated for links to agricultural yield: for example, de la [START_REF] De La Casa | Climate change and its impact on agricultural potential in the central region of Argentina between 1941 and 2010[END_REF] showed a significant correlation between the AMO index and corn and soybean in Argentina. To date, most investigations of the effect of large-scale oceanic/atmospheric variables on agriculture in Europe have been restricted to ENSO (e.g. [START_REF] Cantelaube | Influence of climate variability on European agriculture -analysis of winter wheat production[END_REF][START_REF] Capa-Morocho | Crop yield as a bioclimatic index of El Niño impact in Europe: crop forecast implications[END_REF], 2016a,b;[START_REF] Ceglar | Linking crop yield anomalies to large-scale atmospheric circulation in Europe[END_REF]TNA phenomena (Capa-Morocho et al., 2016b). Considering SSTs and additional observed weather inputs, the Crop Growth and Monitoring System, developed by the European Commission Joint Research Center within the Monitoring Agricultural Resources (JCR-MARS) activities, provides short-term (in-season) yield forecasts of the main food crops in Europe [START_REF] Vossen | Early crop yield assessment of the E.U countries: the system implemented by the Joint Research Center[END_REF][START_REF] Lazar | Methodology of the MARS crop yield forecasting system[END_REF][START_REF] De Wit | Spatial resolution of precipitation and radiation: the effect on regional crop yield forecasts[END_REF][START_REF] Pagani | Improving cereal yield forecasts in Europe -The impact of weather extremes[END_REF]. For grasslands, biomass production is monitored from remote sensing products, but Europe-wide forecasts are not provided yet (MARS, 2019).

Permanent grasslands are agricultural lands used for cultivating "grasses or other herbaceous forage", through either self-seeding or sowing operations, not included in the crop rotation of the farm in at least five years [Reg. EU 2017/2393]. These ecosystems hold a remarkable part of Europe's biodiversity and provide a wide variety of natural resources for human consumption (e.g. livestock, meat and dairy products). Around 20% of the total EU-28 territory is covered by grasslands. Grasslands provide 25% of total food intake for European livestock and play a key role in the greenhouse gas (GHG) budget [START_REF] Hörtnagl | Greenhouse gas fluxes over managed grasslands in Central Europe[END_REF]. In France, permanent grasslands cover 36% of agricultural land area and this number can be as high as ∼60% in some regions like the Massif Central [START_REF] Graux | Ensemble modeling of climate change risks and opportunities for managed grasslands in France[END_REF]. This region, which has the highest number of dairy farms in the country, significantly contributes to keep France ranked amongst the main EU-28 producers of bovine meat and milk (Eurostat, 2019). Considering the importance of the livestock sector in the region, and its dependency on forage resources for animals, there is a great need to develop forage prediction tools to assist farmers in managing risk for grassland-livestock production systems. Establishing a framework for yield prediction matching grassland management, weather and climate is particularly relevant in areas like the Massif Central of France (focus of this study), whose economy is vulnerable to climate variability and change (e.g. [START_REF] Virto | Soil degradation and soil quality in western Europe: current situation and future perspectives[END_REF][START_REF] Hamidov | Impacts of climate change adaptation options on soil functions: a review of European case-studies[END_REF].

The impact of climate on grassland productivity has been evaluated through field campaigns with in-situ meteorological measurements [START_REF] Black | Development and growth characteristics of caucasian and white clover seedlings, compared with perennial ryegrass[END_REF][START_REF] Dürr | Ranges of critical temperature and water potential values for the germination of species worldwide: contribution to a seed trait database[END_REF][START_REF] Chen | The signature of sea surface temperature anomalies on the dynamics of semiarid grassland productivity[END_REF] and remote sensing products [START_REF] Kawabata | Global monitoring of interannual changes in vegetation activities using NDVI and its relationships to temperature and precipitation[END_REF][START_REF] Ding | Temporal interpolation of satellite-derived leaf area index timeseries by introducing spatial-temporal constraints for heterogeneous grasslands[END_REF][START_REF] Liu | Spatiotemporal dynamics of grassland aboveground biomass on the Qinghai-Tibet Plateau based on validated MODIS NDVI[END_REF]. Overall in the mid-latitudes, mild temperatures and generous precipitation are linked to increased biomass production [START_REF] Menzi | Relationship between climatic factors and the dry matter production of swards of different composition at two altitudes[END_REF]. Nevertheless, these factors can notably vary depending on local climatic conditions (e.g. semi-arid, humid, mountainous; [START_REF] Nippert | Intra-annual rainfall variability and grassland productivity: can the past predict the future?[END_REF][START_REF] Yang | Relationship between variability in aboveground net primary production and precipitation in global grasslands[END_REF][START_REF] Carlyle | Response of grassland biomass production to simulated climate change and clipping along an elevation gradient[END_REF][START_REF] Jones | The nitrogen, carbon and greenhouse gas budget of a grazed, cut and fertilised temperate grassland[END_REF]. Most of these studies have provided high-quality grassland data but were limited by the short periods of time over which the analyses were carried out. A remarkable exception can be found in [START_REF] Chen | The signature of sea surface temperature anomalies on the dynamics of semiarid grassland productivity[END_REF], which presented forage yield data over 1939-2016. This is an issue because extended timeseries are scarce but essential to identify robust statistical links between climate variability and yield responses in any agricultural system [START_REF] Craine | Timing of climate variability and grassland productivity[END_REF][START_REF] Capa-Morocho | Crop yield as a bioclimatic index of El Niño impact in Europe: crop forecast implications[END_REF]2016a,b). As an alternative, crop simulation models can be used to generate long-term timeseries of grassland performance data.

Grassland primary production and GHG budgets have been quantified in modeling studies based on observational flux data at specific locations in Europe [START_REF] Ma | Regional-scale analysis of carbon and water cycles on managed grassland systems[END_REF], and world-wide [START_REF] Ehrhardt | Assessing uncertainties in crop and pasture ensemble model simulations of productivity and N2O emissions[END_REF]. Such studies have substantiated the use of grassland biogeochemical models at local, regional and continental scales (Vuichard et al., 2007a) and shown their potential to represent grassland systems under a variety of conditions [START_REF] Brilli | Review and analysis of strengths and weaknesses of agro-ecosystem models for simulating C and N fluxes[END_REF].

Additional opportunities are nowadays offered by high-resolution long-term datasets of observational atmospheric records (e.g., ERA-WATCH, SAFRAN; [START_REF] Weedon | The watch forcing data 1958-2001: a meteorological forcing dataset for land surface-and hydrological models[END_REF][START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the SAFRAN system[END_REF], which provide support to model intercomparison initiatives (e.g., MACSUR, AgMIP, CN-MIP or ISIMIP; [START_REF] Frieler | Assessing the impacts of 1.5 °C global warming -simulation protocol of the inter-sectoral impact model intercomparison project (ISIMIP2b)[END_REF]. It is now possible to feed grassland-specific numerical models with long-term weather observations and perform simulations under potential conditions, i.e. optimal management (N fertilization, number of cuts) and biotic stresses effectively controlled (no pests and diseases). As a result, longterm purely climate-dependent timeseries of potential forage yield can be generated and analyzed with grassland models. This can be done at the cost of downgrading the quality of outputs (i.e. data can be coarse grained for simulation purposes, e.g. [START_REF] Hoffmann | Impact of spatial soil and climate input data aggregation on regional yield simulations[END_REF]. The issue here is to ensure that this downgrading is still adequate for delivering operationally relevant CS.

Based on this novel modeling methodology, a framework was developed in this study to support CS applications and tools with a focus on grassland ecosystems. In particular, we explored the potential links between climate variability, teleconnection patterns and harvested forage in a permanent grassland located in the Massif Central of France. The aim was to provide evidence that skillful seasonal forecasts of forage productivity could be attainable using only oceanic and atmospheric predictors.

Data

Grassland site

The permanent grassland system in this study is located in Laqueuille (45°38′ N, 2°44′ E, 1040 m a.s.l.), in a west-facing leeward slope of the French Massif Central (Fig. 1a andS1a). The grassland plant community is a grass-clover mix dominated by perennial ryegrass (Lolium perenne L.) and white clover (Trifolium repens L.). The soil is a 1 m deep Andosol with 11% carbon and 18% organic matter in the 0.1 m topsoil (they rapidly decline further down) and silt loam texture (26% sand, 20% clay). Soil field capacity, permanent wilting point and saturation are, respectively, 0.36, 0.22 and 0.53 m 3 m -3 [START_REF] Klumpp | Long-term impacts of agricultural practices and climatic variability on carbon storage in a permanent pasture[END_REF]. The grassland is a long-term observation system for research and experimentation managed by INRA since 2003(SOERE, 2019).

Gridded oceanic and atmospheric databases

Gridded oceanic and atmospheric databases were used to analyze climate variability at global, European and country (France) scales. All the information regarding these datasets is summarized in Table S1.

Meteorological observations at the site

An INRA meteorological station at the grassland site provides hourly values of surface air temperature, relative humidity, precipitation, wind speed and incident solar radiation. Due to its mountainous location, the grassland site presents sub-arctic climatic conditions (Dfc Köppen-Geiger classification), with limited sunlight, cold temperatures and abundant precipitation year-round (cf. Table 1 and Fig. 1b). The station records cover the period 1996-2015, but the extensive data gaps in the record of solar radiation only made it possible to perform grassland simulation over eight years (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Instead, we referred to the SAFRAN atmospheric re-analysis (cf. Table S1) to obtain long-term meteorological observations of the area. To this aim, the timeseries of the closest SAFRAN grid cell (central point) to the simulated grassland were selected (3.64 km away, similar altitude). To allow comparison between in-situ and SAFRAN measurements, the maximum overlapping period of both datasets was selected for all weather variables (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), except for surface incident solar radiation (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Table 1 provides the mean and standard deviation (SD) values of both datasets.

On average, a satisfactory agreement was observed between station and SAFRAN data (Table 1). Despite its generally good agreement with Laqueuille station data and its validation throughout France [START_REF] Quintana-Seguí | Analysis of near-surface atmospheric variables: validation of the SAFRAN analysis over France[END_REF], SAFRAN has several systematic biases that deserve further consideration: it slightly overestimates srad, Pr, and T, while slightly underestimating wind speed (wspd). According to Fig. 1b, biases in Pr and wspd are larger during the winter months, and it can be logically assumed that their potential impacts on soil water balance, evapotranspiration and simulated forage yield should be limited, considering that grassland biomass production is minimal from November to March due to low temperatures. For srad, the largest biases arise during spring and summer. These differences could be associated with upslope or local radiation fogs at the grassland site (with plants acting as an important source of water vapor), not adequately captured by SAFRAN. For standard deviation, lower values are obtained in SAFRAN for T, relative humidity and wspd at daily and interannual timescales compared to the INRA station (Table 1). Contrastingly, srad variance is overestimated in SAFRAN. Indeed, in the computation of anomalies for this study, all means were removed, thus minimizing the potential influence of biased data in the final results. As a consequence, only the ), maximum and minimum daily temperatures (red/blue -°C), daily precipitation (cyan -mm day -1 ) and wind speed (magenta -m s -1 ) from the in situ meteorological station (solid lines and bars) and SAFRAN closest grid cell data to Laqueuille (dashed lines and stems). Period 1996-2015 except for solar radiation (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).

SAFRAN weather data (57 years) were utilized as an input for grassland modeling in this study, discarding the weather station data owing to their short temporal extent (8 years).

Methods

Statistical methods and teleconnection indices

Climate anomalies were calculated by subtracting the mean value from a timeseries of weather data and then dividing the results by its standard deviation (standardization process). For seasonal/monthly anomalies, the mean of the corresponding season/month was thus subtracted to remove the impact of the seasonal cycle.

To separate contributions from interannual and lower frequency variability (e.g., multi-decadal), a zero-phase spectral 10th-order Butterworth filter with 11-year cut-off period was utilized in high-pass and low-pass modes. The choice of 11-year cut-off period was based on previous studies with a similar focus on interannual and lower-frequency climate variability [START_REF] Gómara | Abrupt transitions in the NAO control of explosive North Atlantic cyclone development[END_REF]. The Butterworth method [START_REF] Rabiner | Theory and Application of Digital Signal Processing[END_REF] has the advantage that it efficiently removes trends in high-pass mode and minimizes undesired edge effects compared to analogous spectral filters (e.g., Lanczos; cf. Supplementary Fig. S1b for an example of surface radiation timeseries).

The methodologies applied to calculate the teleconnection indices are specified in Table S2. A regression of interannual mslp anomalies on the NAO index is provided in Fig. 2a (NAO+ pattern -shadings). SST patterns associated with the Niño3.4, TNA, MED and AMO indices are provided in Fig. 2b (positive phases). Please note that they appear plotted all together in Fig. 2b just for space saving and not to indicate any potential coincidence in time.

Student-t tests that account for the autocorrelation of the timeseries through the calculation of effective degrees of freedom were used for hypothesis testing [START_REF] Bretherton | Effective number of degrees of freedom of a spatial field[END_REF].

The pasture simulation model (PaSim)

Grassland simulations were performed with the Pasture Simulation model (PaSim). PaSim is a deterministic, biogeochemical grassland-plot model (initially developed by [START_REF] Riedo | A pasture simulation model for dry matter production and fluxes of carbon, nitrogen, water and energy[END_REF] incorporating climate data, soil properties, vegetation characteristics, livestock (dairy/beef cattle and sheep) and management, and operating at a point in space on an hourly time step. PaSim simulates multi-year growth of perennial species. It considers the vegetation cover as a single-plant community (although a fixed percentage of legumes can be set to simulate symbiotic N fixation). The model is suited to assess the effect of management practices in terms of fertilization, grazing intensity and duration, and cutting frequency.

The model calculates water, C and N pools and fluxes. Photosynthetic-assimilated C is allocated at each time step to four pools (roots, leaves, stems/sheaths and ears). For each organ, the biosynthesis pathway implies a transition on four age classes from newly produced tissue until senescence. Animal milk production, enteric methane emissions, nitrogen emissions from nitrogen input flows (including returns to soil from animal excreta), and ecosystem respiration are outflow fluxes. Accumulated above ground biomass, if not mown or grazed, enters litter reservoirs. The litter is evenly distributed into the whole soil profile, segregated into its structural and substrate components. The soil organic matter also differentiates between active, slow and passive pools with different decomposition rates according to firstorder kinetics. N inputs to the soil include atmospheric deposition, biological fixation by legumes and fertilizer addition. Losses of N occur via pathways that include nitrate leaching, ammonia volatilization and Table 1 Daily mean values of meteorological variables from the INRA in situ station and SAFRAN closest grid cell to Laqueuille. For clarity, the total annual values of precipitation are incorporated. Standard deviation values of all variables are also provided. All values are obtained from the 1996-2015 period except surface solar radiation (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Missing values from station data represent at most 0.006% of total for the periods and variables selected. gaseous emissions through microbial conversion of ammonium and nitrate (such as nitrification and denitrification). In the ground process scheme, the soil profile is divided into six layers. In each layer, the kinetics for N flows are driven by C/N ratios with the following values; structural: 150; metabolic: 10; active: simulated in the range 3-14; slow: simulated in the range 12-20; passive: simulated in the range 6-8. Soil temperature and moisture are simulated as a function of soil physical properties and plant water use in each layer of soil.

The model is enabled for gridded simulations [START_REF] Vital | High-performance computing for climate change impact studies with the pasture simulation model[END_REF][START_REF] Eza | An open platform to assess vulnerabilities to climate change: an application to agricultural systems[END_REF] over large-scale areas and for climate change-impact studies [START_REF] Graux | Development of the pasture simulation model for assessing livestock production under climate change[END_REF][START_REF] Graux | Ensemble modeling of climate change risks and opportunities for managed grasslands in France[END_REF][START_REF] Lardy | Numerical methods for project vulnerability analyses with guidelines for use[END_REF][START_REF] Ma | Regional-scale analysis of carbon and water cycles on managed grassland systems[END_REF][START_REF] Pulina | Modeling pasture production and soil temperature, water and carbon fluxes in Mediterranean grassland systems with the pasture simulation model[END_REF]. PaSim is also part of multi-modeling exercises in the frame of international initiatives [START_REF] Sándor | Uncertainty in simulating biomass yield and carbon -water fluxes from grasslands under climate change[END_REF](Sándor et al., , 2016a(Sándor et al., ,b, 2017(Sándor et al., , 2018a;;[START_REF] Ehrhardt | Assessing uncertainties in crop and pasture ensemble model simulations of productivity and N2O emissions[END_REF]. Sándor et al. (2018b) improved the model to account explicitly for thermal acclimation of plants. In addition to SAFRAN weather data, NOAA Mauna Loa CO 2 annual mean data (NOAA-ESRL, 2019) were used as input for PaSim.

Experimental setup: grassland management and model calibration

Alternative grassland management options for modeling were designed with different levels of complexity. Specific details of all performed simulations are provided in Table 2 and Fig. 3.

All runs simulated potential forage yield, whose definition closely mirrors the one used for arable crops, i.e. "the yield of a cultivar when grown in environments to which it is adapted, with nutrients and water non-limiting and with pests, diseases, weeds, lodging, and other stresses effectively controlled" [START_REF] Evans | Yield potential: its definition, measurement, and significance[END_REF]. In this way, the gap between observed and simulated yield is sensitive to the optimal management conditions (e.g. fertilization, irrigation, cuts), among other potential factors (e.g. model errors).

The choice to use purely simulated forage data rather than empirical data was due to: (i) the short time period covered by the available observational data (∼15 years), which hampered inferring any robust conclusion in the context of this work, as generally periods of at least 30 years are desirable to establish robust links and analyze causality between year-to-year climatic and non-climatic variables [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF]IPCC, 2014); and (ii) observed data were influenced by other factors like pests, technical development or research priorities, which obscured the purely climatic origin of the vegetation response that this study aimed to identify. Simulated data has helped to overcome these issues in previous studies, as the most convenient to reveal such interactions [START_REF] Capa-Morocho | Crop yield as a bioclimatic index of El Niño impact in Europe: crop forecast implications[END_REF], 2016a,b).

The choice of principally simulating mown grasslands in an environment where half of the land is grazed [START_REF] Rapey | Multiple forms of mountain-area dairy farms converting from milk to meat[END_REF][START_REF] Estel | Combining satellite data and agricultural statistics to map grassland management intensity in Europe[END_REF] instead of using actual empirical grazing data also responded to the need to generate purely climate-dependent timeseries of forage production. For this purpose, the simplest management conditions (cf. MSB in Table 2 and Fig. 3a, which does not include any cutting event) were first evaluated.

Subsequently, management complexity was increased via more frequent cutting events (e.g., 2C, 3C; Fig. 3b andc). Finally, we considered the most complex scenario, denoted 'optimal management' or OM (Fig. 3d), which mimics mowing grassland management generally used by farmers in Europe [START_REF] Vuichard | Estimating the greenhouse gas fluxes of European grasslands with a process-based model: 2. simulations at the continental level[END_REF]. Mown grasslands are not uncommon in the French Massif Central area, where conserved forage (silage and hay) is often utilized by dairy farms to ensure cattle feeding requirements in winter (e.g. [START_REF] Violleau | L'ensilage d'herbe dans les systèmes bovins laitiers du nord du Massif Central (cas du Puy-de-Dôme)[END_REF][START_REF] Baumont | Transformation des plantes au cours de leur conservation et conséquences sur leur valeur pour les ruminants[END_REF].

Regarding the grassland model setup, a regional 'European' calibration, obtained by calibrating the model against multiple locations (either grazed and/or mown) in Europe [START_REF] Ma | Regional-scale analysis of carbon and water cycles on managed grassland systems[END_REF], was applied for the site of Laqueuille. The fact that Laqueuille is a grazed system and not mown as in our simulation, enforced the use of regionally-adjusted parameters instead of site-specific calibrated ones, thus reducing simulation accuracy in a way that could be acceptable to ensure a general for each application (higher N amounts provide same results). Irrigation disabled (as in the rest of simulations if not stated otherwise). Automatic cuts disabled. Although no forage harvests take place in this simulation, the yield is assumed to be equal to the maximum shoot biomass attained in the middle of the growing season. As an example, Fig. 3a provides the simulated potential shoot biomass of year 1982 (maximum around day 175 -24 June).

6.87 0.9 2 Cuts

2C

Two annual forage cuts. The first when the value of maximum shoot biomass is reached (same as in MSB, with varying cutting dates between years). The second when a second maximum after regrowth is attained. Potential forage yield value of each year is retained (Fig. 3b green line). Additional to the N applications on days 14 and 120 (160 kg N ha -1 each), a third N application is triggered after the first annual cut.

OM

The optimal management routine from [START_REF] Vuichard | Estimating the greenhouse gas fluxes of European grasslands with a process-based model: 2. simulations at the continental level[END_REF] is utilized to simulate the potential forage yield under the possible closest circumstances to reality. Cutting dates vary each year and the annual potential forage yield values are retained. The first cut of each year is triggered right before plant senescence. The subsequent cuts take place whenever plant growth rate declines during 10 consecutive days and there is enough available forage to collect (> 1 t DM ha -1

) after at least 30 days of regrowth.

Fertilization events are optimized based on a N nutrition index (N application amounts reduced to minimum non-limiting forage growth levels). In Fig. 3d, the management for year 1982 is also provided as an example.

7.60

1.5 applicability of the modeling framework. By accepting the extrapolation principle from [START_REF] Ma | Regional-scale analysis of carbon and water cycles on managed grassland systems[END_REF], we held it valuable to extrapolate towards a mown grassland in Laqueuille, where rich site-specific soil, plant and weather data were also available to back the study. The choice of Laqueuille for this pilot study was also appealing because this site is conveniently located in the main grassland region of France, the Massif Central.

Statistical seasonal forecast model of potential forage yield

A statistical seasonal forecast model of potential forage yield at Laqueuille was developed. A stepwise approximation was followed for adding or removing terms in the multi-linear regression model. The method automatically searches for sets of oceanic/atmospheric predictive terms improving forecast skill through forward selection and background elimination. At each step, the p-value of an F-statistic is utilized to select 'in' and 'out' predictors [START_REF] Draper | Applied Regression Analysis[END_REF]. Next, a leave-one-out cross-validation method was applied, where the model was trained from predictor/predictand data combinations from all years except the year being predicted at each training stage [START_REF] Allen | Mean square error of prediction as a criterion for selecting variables[END_REF].

Results and discussion

This section first illustrates the temporal evolution of grassland productivity under alternative management options and analyzes possible correlations between forage production and biophysical factors (weather variables and CO 2 concentration; Section 4.1). Then, the contributions of trends/low-frequency (Section 4.2) and high-frequency (Section 4.3) climate variability on forage yield are presented separately. Results and discussion of the developed statistical seasonal forecast model are provided in Section 4.4.

Temporal evolution of grassland productivity

With higher number of cuts (with or without irrigation; 3CI/3C), simulations present higher mean productivity than with limited cuts (2C) or maximum shoot biomass annual values (MSB; Fig. 4a andTable 2).

However, the optimal management simulation (OM) does not generate the highest potential forage yield, owing to internal biophysical constraints on the simulated forage production (e.g., water soil bearing capacity may not adequately support loads, or cutting events after 30 days of regrowth are only triggered whenever sufficient forage yield is available for harvesting, > 1 t DM ha -1 ).

The standard deviation values (Table 2) appear to arise from different temporally-dependent contributions (Fig. 4a). The first and most evident contribution is associated with trends and/or potential low frequency oscillations in forage yield. For instance, MSB increases by 29% (+1.74 t DM ha -1 ) between 1959-1979 and 1995-2015. The second contribution is due to interannual yield changes, which is more evident under more intensive management options (3C, OM). To analyze possible climate-related causes for these changes, Fig. 4b provides the yearly averaged standardized anomalies of observed climatic conditions at the grassland together with yield for the MSB simulation, which is the simplest management scheme. Anomalies of surface radiation (srad), maximum and minimum temperatures (Tmax/ Tmin) and CO 2 concentration appear to covary with yield, both in terms of year-to-year and longer-term changes.

The correlation values between these climatic variables (plus precipitation) and yield are provided in Table 3 (left values), srad being the main driver of forage yield variability in all simulations. For completeness, correlation coefficients were also calculated between yield and other climatic variables also utilized by PaSim to estimate water balance (water vapor pressure and wind speed), but they did not provide anything worth commenting on.

These results are consistent with the Alpine-like climate of the grassland site, where incident solar radiation is limited due to frequent rainfall (cf. Table 1 and Fig. 1b) and temperatures are low due to the high altitude (1040 m a.s.l.). Higher values of CO 2 , srad and Tmax/ Tmin (all highly correlated with yield) are known to stimulate plant photosynthesis and promote shoot biomass growth in Alpine environments [START_REF] Menzi | Relationship between climatic factors and the dry matter production of swards of different composition at two altitudes[END_REF][START_REF] Ainsworth | The response of photosynthesis and stomatal conductance to rising [CO 2 ]: mechanisms and environmental interactions[END_REF]. Logically, applying frequent forage cuts during the year (3C, OM) introduces higher variability in simulated grassland productivity than just measuring the maximum shoot biomass attained in the middle of the growing season of each year (standard deviation values in Table 2), which typically occurs in June (MSB; cf. Fig. 3a). This becomes especially evident after the decade of the 1990s, when atmospheric and climatic conditions considerably improved (higher CO 2 , Tmax/Tmin

Table 3

Correlation values between annual averaged anomalies of potential maximum shoot biomass / forage yield and SAFRAN climatic conditions at the grassland (period 1959-2015): incident solar radiation (srad), maximum/minimum daily temperature (Tmax/Tmin), precipitation (Pr) and CO 2 concentration. Values to the left/ center/right are obtained from raw/low-pass/high-pass filtered timeseries. Boldface indicates 95% confidence interval (t-test). and srad; Fig. 4b), raising the yield ceiling and making strategic management a much more important driver of productivity (Fig. 4a).

For precipitation, the correlation with yield is, in general, slightly negative and not statistically significant (p > 0.05). Again, this is consistent with the high rainfall and low radiation at the site. First, rainfall abundance prevents relevant water stress on grassland growth. Second, yield increases when radiation and surface temperature increase, which in turn are usually associated to lower cloudiness and precipitation. Both interactions together explain the precipitation-yield response.

Impact of trends and low-frequency climate variability on grassland productivity

The correlations between low-pass filtered timeseries of yield and annual climatic anomalies are provided in Table 3 (central values). Positive and statistically significant (p < 0.05) correlations can be observed between yield and srad, Tmax/Tmin and CO 2 , and negative for Pr. This is particularly evident between srad and the different biomass simulations (0.95-0.98 correlation values). Thus, the long-term changes of the variables explaining a higher fraction of yield variability, srad, Tmax and CO 2 , were subsequently analyzed.

In Fig. 4c, the raw timeseries of MSB yield is provided together with low-pass filtered srad and Tmax. Apart from the filtering, the main difference between Fig. 4c andb is that srad and Tmax timeseries (in gray/red) are based on spatial averages over mainland France [5°W-8°E , 42°N-52°N] and not just Laqueuille. This indicates that the increase in srad and Tmax between 1959-1979[START_REF] Vossen | Early crop yield assessment of the E.U countries: the system implemented by the Joint Research Center[END_REF]-2015 (Fig. 4c (Fig. 4c) was associated to a large-scale feature.

Regarding the potential causes for the large-scale increase in srad over France, several authors have pointed to the global 'dimming and brightening' effect related to anthropogenic emission of aerosols [START_REF] Wild | From dimming to brightening: decadal changes in solar radiation at Earth's surface[END_REF][START_REF] Sánchez-Lorenzo | Dimming/brightening over the Iberian Peninsula: trends in sunshine duration and cloud cover and their relations with atmospheric circulation[END_REF]. However, the level of detail of this study does not allow disclosing the radiative effect of aerosols (absorption or scattering). Instead, multi-decadal changes in cloud cover fraction seem to be the main driver of srad changes (cf. Fig. 4c -gray and magenta lines). The multi-decadal NAO index is overlaid in Fig. 4c because NAO is one of the main drivers of multidecadal cloudiness variability over the North Atlantic [START_REF] Trigo | The North Atlantic Oscillation influence on Europe: climate impacts and associated physical mechanisms[END_REF], due to its influence on low pressure systems' trajectories [START_REF] Ulbrich | A shift of the NAO and increasing storm track activity over Europe due to anthropogenic greenhouse gas forcing[END_REF][START_REF] Gómara | On the relation between explosive cyclones affecting Europe and the North Atlantic Oscillation[END_REF][START_REF] Gómara | Abrupt transitions in the NAO control of explosive North Atlantic cyclone development[END_REF]. The index may depict a remarkable variability [START_REF] Wang | Decadal variability of the NAO: introducing an augmented NAO index[END_REF][START_REF] Pinto | Past and recent changes in the North Atlantic Oscillation[END_REF], but the observed behavior is not consistent with srad changes. Another teleconnection pattern with significant influence on European climate at these timescales is the AMO (Fig. 2b; [START_REF] Knight | A signature of persistent natural thermohaline circulation cycles in observed climate[END_REF]Wollings et al., 2012). The AMO index (also overlaid in Fig. 4c) does appear to covary with srad and cloud cover changes. Nevertheless, it is still only a brief snapshot of a long-term cycle. A negative AMOcloud cover relation in the North Atlantic, mediated by SST meridional gradients and storm track activity, has been recently pointed out by [START_REF] Vaideanu | Atlantic multidecadal oscillation footprint on global high cloud cover[END_REF].

The positive trends in Tmax/Tmin and CO 2 concentrations (Fig. 4b andc) are in line with global warming trends (IPCC, 2014). Trends of Tmax values show a similar increase at Laqueuille and over mainland France when disaggregated among seasons (Figs. S1c-d). The same results were obtained for Tmin (data not shown).

Investigating further these aspects, though interesting, would deviate too much from the initial objectives of this study, and could be arranged later (e.g., through climate-grassland sensitivity model experiments) as a natural evolution of what is presented here.

Impact of interannual and seasonal climate variability on grassland productivity 4.3.1. Relations between grassland shoot growth rate, forage yield and climate anomalies

To analyze the influence of interannual climate variability on grassland productivity, all considered timeseries in this section were high-pass filtered (if not stated otherwise). The correlation values between climate and forage yield anomalies are shown for all simulations in Table 3 (right values). At interannual timescales, the main drivers of forage productivity are again srad (MSB/3C/3CI simulations) and Tmax/Tmin (3C/3CI/OM). Despite still being statistically significant, the observed correlations are clearly lower at interannual timescales than for the raw timeseries. In addition, correlation coefficients vary among simulations (e.g., srad/Tmax/Tmin values for MSB and 3C/3CI; Table 3).

To identify the seasons of the year in which climate conditions present a stronger influence on forage growth, Fig. 5 (left column) provides the correlation values between monthly climate and shoot growth rate anomalies (MSB, 3C and OM simulations). Results for 2C and 3CI are also available in supplementary Fig. S2. As expected, Tmax/Tmin and srad positive anomalies at the grassland are robustly linked (p < 0.05) with enhanced forage growth rates for long periods of the year. For clarity, the monthly mean shoot growth rates appear overlaid in Fig. 5 (dotted black line in panels). Again, the negative correlation with precipitation is a consequence of the high and homogeneous rainfall, and the positive effect of sunnier and warmer weather conditions on forage growth at this mountainous location. In addition, the periods in which these relations are stronger are sensitive to the timing of cutting events and seasonal changes in climate conditions. For MSB (Fig. 5a), the period of the strongest climate-growth rate links spans from January to July. This is because maximum shoot biomass values tend to occur between June and August. For 3C (Fig. 5c), as the three annual cuts span from June to November, the influence of Tmax/ Tmin and srad on growth rate is extended until the end of the year (from January to November).

During the summer months this influence appears slightly reduced though. This behavior is also observed in the 2C simulation (Fig. S2a) and might be consistent with increased water demand during summer, thus requiring a more balanced equilibrium between Tmax/Tmin/srad and Pr effects. This hypothesis is confirmed in Fig. S2c, where the 3CI simulation (which performs optimal irrigation removing Pr dependency) shows higher correlation values with srad and T during the summer months compared to 3C (Fig. 5c). For OM, results are very similar to the 3C simulation (Fig. 5c ande).

Next, the anomalies of total annual forage yield (1 value per year) were correlated against monthly climatic anomalies (12 values per year). Therefore, for each simulation, the same yield timeseries was correlated against the different monthly climatic anomalies (Fig. 5, right column). In general, the links between forage yield and climate anomalies are notably weaker compared to forage growth rates. This is because correlation values are no longer based on concurrent monthly anomalies. Instead, yield values provide combined information of climate conditions over long periods of the year and the efficiency of grassland management to materialize increased forage growth rates into actual yield. Thus, simulations indicate that the most intensively managed grasslands are better linking srad/Tmax/Tmin and yield anomalies, especially during spring and fall. This can be clearly observed when comparing 3C/OM/3CI simulations (Figs. 5d/f/S2d) with MSB/2C ones (Figs. 5b/S2b). During the summer months similar results point to increased Pr dependency, caused by increased water demand habitual of this season (e.g. 3C; Fig. 5d). This increased Pr dependency is logically more evident for correlations based on longer temporal windows (compare blue line in Fig. 5c -monthly-versus d -annual).

Relations between forage yield and large-scale monthly climate anomalies

To confirm whether the climate conditions fostering grassland productivity were specific to the site location (e.g., absence of upslope fogs) or related to any large-scale atmospheric circulation pattern, the correlation between forage yield at Laqueuille and SAFRAN climate anomalies was assessed over mainland France (Figs. 6 and 7;left and central columns). The choice of variables (Tmax/Tmin/srad), simulations (3C/OM) and months (January/February/April/September) responds to Fig. 5 outcomes: 3C and OM simulations provide the strongest links between meteorological surface variables and yield during the selected months. Besides, 3C and OM simulations are respectively the most productive and closest to real-life ones (Table 2). As it can be observed, statistically significant correlations (in stippling) extend over vast areas. This implies that warmer temperatures and enhanced srad that produce higher forage yields at Laqueuille are generally present at the same time over large areas of mainland France. This is particularly more evident for 3C, as the regions of significant positive correlations are larger than for OM, for the latter being, in some cases, even of opposite sign, but not statistically significant (Fig. 7e). For a better observed mslp structure as the spatial correlation between the composite of mslp maps from Figs. 6 and 7 (eight members in total; pattern in Fig. 2a -contours) and mslp high-pass filtered monthly anomalies for the period 1959-2015. Hereafter, we refer to it as the mslp-index (cf. also Table S2).

The correlation value between NAO and mspl-index is 0.81 (p < 0.01), while the correlation with the East Atlantic pattern, a similar mode to NAO but with centers of action shifted south, is 0.25 (p < 0.1) [START_REF] Barnston | Classification, seasonality and persistence of lowfrequency atmospheric circulation patterns[END_REF]their Fig. 9d). These results indicate that the identified mslp pattern appears to enhance Tmax/ Tmin and srad over large areas of France, including the Massif Central, which in turn promotes forage yield at Laqueuille. This is consistent with positive mslp (anticyclonic) anomalies over the area, which are typically associated with more stable weather; sunnier and warmer conditions during daytime. For completeness, the same maps covering the North Atlantic area are provided for all months for 3C and OM simulations in Figs. S5 andS6, respectively.

Regarding SSTs, Figs. 6 and 7 also depict a positive robust connection between enhanced 3C/OM yield and a warming of the marine waters near France (Mediterranean, North Sea and far eastern North Atlantic). This is especially evident in the months of February and September (Figs. 6 and 7f,l). As described in the introduction, these SST conditions are typically linked to higher temperatures over Europe [START_REF] Feudale | Role of Mediterranean SST in enhancing the European heat wave of summer 2003[END_REF][START_REF] García-Serrano | Large-scale atmospheric response to eastern Mediterranean summer-autumn SST anomalies and the associated regional impact[END_REF].

To test the sensitivity of results to the choice of atmospheric reanalysis (NCEP, ERA-40, ERA-Interim; cf. Table S1) and provide a global view on the relation between mslp/SST and 3C/OM yield anomalies, in Figs. S7-S8 the analogous correlation maps are provided. In all cases, the 3C/OM yield timeseries are the same (cf. Table 2 forced by SAFRAN) and the only factors that vary are the mslp anomalies from each re-analysis database and the period utilized to compute the correlations (cf. Table S1). For simplicity, in Fig. 8 the main outcomes from Figs. S7-S8 panels are outlined.

A first conclusion is that correlation patterns of mslp are quite similar among re-analyses. A second conclusion is that additional SST patterns to MED also appear to influence simulated forage yield, but only for the 3C simulation. This is the case between 3C yield and TNA SSTs in January (Fig. 8a-c) and with pacific equatorial SSTs (ENSO) in May (Fig. 8g-i). More precisely, supplementary Figs. S7 and S9 reveal that ENSO-3C yield connection takes place between April and July (Fig. S7j-u) and that TNA influences yield from September of the previous year (hereafter y-1; Fig. S9i-l) to January of harvest year (y0; Fig. S7a). Additionally, these SST-yield relationships appear to evolve in time. For instance, the influence of ENSO on 3C yield is much stronger in ERA-Interim (period 1979-2015) than in NCEP (period 1959-2015; compare Fig. 8g andi).

The differences observed between 3C (Fig. S7) and OM (Fig. S8) simulations may be caused by the much more complex management and limitations (soil bearing capacity, sufficient yield for harvesting) in the OM simulation. The fact that 3C cuts are fixed in time every year and OM ones not could also contribute to sharpen ENSO/TNA SST climatic impacts on yield. For instance, if ENSO influenced climate conditions during summer, this signal could be captured by the 2nd cut in 3C (always covers the period June to August). For OM, as cuts are variable in time, this relation could get attenuated.

As a summary of the relations identified so far on NCEP (the longest re-analysis), in Fig. 9 the lagged correlation values of yield versus mslp (NAO, mspl-index) and SST (TNA, ENSO, MED) teleconnection indices are provided for 3C (Fig. 9a) and OM (Fig. 9b). As expected, the msplindex is significantly correlated with 3C yield on February, April and September y0 (p < 0.05). Similar results are obtained for NAO and MED on the same months. Niño3.4 is significantly correlated with yield in May and June y0 and TNA from September y-1 to January y0. For OM, only the links between yield and mspl-index/MED remain for February and September.

Influence of teleconnection patterns on forage yield variability

To assess the individual influence of the identified teleconnections on the observed climatic and yield anomalies at Laqueuille (Figs. 5678) and propose a bio-physical mechanism for forage growth, a reverse methodology was applied.

For this purpose, the indices of the teleconnection patterns were directly correlated with monthly anomalies of climate variables at Laqueuille during the previous and following months over France (SAFRAN) and the whole North Atlantic (NCEP). The selected teleconnections for analysis, based on Fig. 9a,b results, were: (i) mspl-index in February and September (3C and OM); (ii) MED in February and September (3C and OM); (iii) ENSO in May-June (3C); and (iv) TNA from September y-1 to January y0 (3C).

Results for (i) are provided in Fig. 9c,d. In Fig. 9c correlations between February mspl-index and January to June climate anomalies indicate that this teleconnection pattern is associated with increased Tmin/Tmax and srad over Laqueuille from January to May. In the same figure the correlations of Feb/mspl-index with the NAO and itself for the previous and following months (black solid/dotted lines) indicate that the persistence of this atmospheric pattern is of just 2 months (February and March).

For completeness, the correlation maps between Feb/mspl-index and concurrent monthly anomalies of mslp (NCEP; red/blue contours), SST (HadSST; shadings), moisture flux at 850 hPa (NCEP; arrows), Tmax (SAFRAN; shadings over France) and srad (SAFRAN; black contours over France) are shown in Fig. 9d. For the calculation of moisture flux, the zonal and meridional wind anomalies at 850 hPa were multiplied by the same pressure level specific humidity values (in m g s -1 kg -1 ). As expected, the anticyclonic conditions associated with msplindex over central France promote higher temperatures and srad (Fig. 9d) and seem responsible for the enhanced forage growth rates at Laqueuille between February and May (Fig. 9c -green line).

In Fig. 9e,f the analogous results are provided for the February MED index. In this case, the autocorrelation of the Feb/MED SST index for the previous and following months indicates a higher temporal persistence of the anomalies in time (magenta line in Fig. 9e). This is due to the high thermal capacity of SSTs [START_REF] Bjerknes | Atmospheric teleconnections from the equatorial pacific[END_REF]. A warmer Mediterranean is also associated with increased temperatures and precipitation at Laqueuille in January and February and higher forage shoot growth rates from January to April (Fig. 9e).

The associated large-scale correlation maps with Feb/MED (Fig. 9f) indicate that a warmer Mediterranean seems to be associated with enhanced temperatures over central Europe and a dipolar mslp structure over the North Atlantic. This mslp structure resembles the East Atlantic pattern and is suggested to be a potential driver of interannual SST winter anomalies in the Mediterranean [START_REF] Skliris | Long-term sea surface temperature variability in the Aegean Sea[END_REF][START_REF] Zveryaev | Interannual variability of Mediterranean evaporation and its relation to regional climate[END_REF].

Attending to the timing and persistence of mspl-index and MED in Fig. 9c ande, results suggest that the MED SST anomalies grow earlier (January) but are not necessarily responsible of the mspl-index pattern onset in February (black line in Fig. 9e). However, whenever msplindex is established in February, the associated anticyclonic conditions over the Mediterranean could as well contribute to the persistence of the positive SSTs over the area (positive feedback -magenta line in Fig. 9c). The combined effect of mspl-index and MED patterns seems responsible for the enhanced forage growth rate values in Fig. 9c until May (green line), a behavior which cannot be explained by the effect of the mspl-index alone.

Regarding precipitation, the enhanced values in January/February associated to a warmer Mediterranean (Fig. 9e) could be explained by the increased moisture flux values over France in Fig. 9f (arrows). The positive SST anomalies also present near the European Atlantic coast may enhance surface evaporation and humidity, which is subsequently transported towards France by the anomalous wind flow.

For completeness, the analogous of Fig. 9c-f but for the month of September are provided in Fig. S10 (June to November period). The results for September are highly similar to those of February and indicate that mspl-index and MED teleconnection patterns are as well responsible of increased shoot biomass productivity at Laqueuille in August and September.

The lead/lag correlations for ENSO and TNA indices did not provide clear and persistent statistically significant anomalies of Tmin, Tmax, srad or shoot growth rates during the harvest year on NCEP (not shown). Therefore, a direct joint climate and eco-physiological mechanism for forage growth was hard to establish. Typically, a positive TNA in summer and autumn is associated with NAO-during the next winter and enhanced precipitation over southwestern Europe [START_REF] Okumura | Tropical Atlantic air-sea interaction and its influence on the NAO[END_REF][START_REF] Rodríguez-Fonseca | Evaluation of the North Atlantic SST forcing on the European and North African winter climate[END_REF]. In this line, one potential mechanism for enhanced yield may be associated with increased soil moisture storage during the previous autumn and winter to forage collection, thus ensuring water supply within a certain range for plant growth in the following months [START_REF] Craine | Timing of climate variability and grassland productivity[END_REF]. For ENSO, as the highest correlation values were found in May and June (Fig. 9a), one possible mechanism could be that the associated signal influenced the precipitation regime in Laqueuille during summer and fall [START_REF] Shaman | The seasonal effects of ENSO on european precipitation: observational analysis[END_REF]his Fig. 1c). In this line, several studies have shown that increased precipitation over southern/central Europe occurs under El Niño conditions [START_REF] Mariotti | Euro-Mediterranean rainfall and ENSO-a seasonally varying relationship[END_REF][START_REF] Park | Remote ENSO influence on Mediterranean sky conditions during late summer and autumn: evidence for a slowly evolving atmospheric bridge[END_REF][START_REF] Shaman | The seasonal effects of ENSO on european precipitation: observational analysis[END_REF], enhancing crop productivity over the Iberian Peninsula [START_REF] Capa-Morocho | Crop yield as a bioclimatic index of El Niño impact in Europe: crop forecast implications[END_REF]2016a,b,c). Consequently, water stress conditions during the summer months (detected in Fig. 5c,d) may be minimized and forage growth could be promoted.

Regarding these mechanisms, a much deeper analysis may be needed to establish the physical links between the SST precursors and climate anomalies at Laqueuille. Instead, in the next section an effort was made to forecast 3C and OM grassland annual forage yield values at Laqueuille based on climatic predictors identified so far.

Statistical seasonal forecast of forage yield

Model performance

A statistical multi-regression model was built to predict annual values of yield anomalies for the 3C and OM simulations. The purpose of the model is to forecast, for each year, whether the forage productivity at the grassland will be 'good' (above average) or 'bad' (below average) several months ahead. To this aim, local climatic monthly anomalies and large-scale monthly teleconnection indices from y0 and y-1 were considered as predictors. As in the previous sections, all considered timeseries were high-pass filtered. The considered predictors for the model were: (i) local monthly anomalies at Laqueuille: Tmax, Tmin, srad and Pr; and (ii) large-scale monthly teleconnection indices: NAO, mspl-index, MED, TNA and Niño3.4. The initial choice of predictors was based on Figs. 5-8 outcomes. Next, a stepwise approximation was followed for adding or removing terms in the multi-linear regression model [START_REF] Draper | Applied Regression Analysis[END_REF].

The construction of a prediction model was first carried out for the 3C simulation. As annual cuts of forage are always triggered on days 150 (30 May), 225 (13 August) and 300 (27 October), atmospheric and oceanic predictors were only considered until the month of May y0. Therefore, the forecast windows with cuts are 0 (cut1), 2.5 (cut2) and 5 (cut3) months, with the total annual yield value always measured after cut3. In this sense, the contributions of the three cuts into the total yield annual value are, respectively, 36, 43 and 21% (1959-2015 averages). These numbers indicate a high degree of actual prediction (64% of total forage; sum of 2nd and 3rd cuts) under the current model configuration. Based on the selected variables, the stepwise method returned 6 predictors providing model skill: [TNA September y-1, mspl-index January y0, mspl-index February y0, MED February y0, Tmax Laqueuille April y0 and Niño3.4 May y0]. Next, a leave-one-out cross-validation method was applied.

Fig. 10a provides the 'observed' (simulated by PaSim) and predicted timeseries of forage yield anomalies based on the generated statistical model. The model performs quite well, providing a skill in correlation of 0.55 (p < 0.01, t-test), a root mean square error percentage (% RMSE) of 14.6% and a 71.43% hit rate in determining whether the yield of a given year will be above or below average. Please note that for the latter the baseline prediction skill is 50%. In the form of a tercileforecast (33% superior, central or inferior values), the model returns a 48.21% hit rate, the baseline being 33%. The same analysis was subsequently performed for the OM simulation. In this case, the number of selected predictors by the stepwise method was notably lower: [MED February y0 and Tmin Laqueuille July y0]. This was somewhat expected due to the low number of statistically significant correlation values between OM yield, climatic anomalies (Fig. 5f) and teleconnection indices (Fig. 9b). To improve the balance between the actual forecast window and model performance, the period for potential predictors was extended until July y0. For the OM simulation, cutting dates are variable among years, triggered from April to November. That provides a 0 to 4-month forecast window between atmospheric/oceanic predictors and yield, but also some temporal overlapping between April and July. Indeed, the model could have been improved by introducing yield values from the first cuts of the year as predictors, but this would have gone in essence against the initial objectives of this paper. The results for the regression model after the leave-one-out cross-validation are provided in Fig. 10b. In this case, the skill in correlation is 0.41 (p < 0.01), %RMSE 17.6% and the hit rates of above/below average annual harvest and tercile-forecasts are 67.86% and 44.64%, respectively (50% and 33% baselines).

Analysis of stationarity in model predictability

This section analyzes the evolution in time of the links between forage yield and climate predictors. This was motivated by the results obtained in Fig. 8g-i, where the link between ENSO-3C yield appeared to be sensitive to the period considered. For this purpose, the 21-year centered running correlations between the 'observed' and 'predicted' 3C/OM yield timeseries are also overlaid in Fig. 10a,b (blue line; p < 0.05 in circles).

For 3C (Fig. 10a), a much higher skill in correlation was found between 1980 and 2015 (max. 0.75) compared to 1959-1980 (min. 0.1). To analyze the causes underlying this behavior, in Fig. 10c the same running correlations were applied between 3C yield and model predictors separately. Results indicate that the main source of 3C yield predictability seems to be ENSO and that this link is not stationary in time. During the periods in which the link between ENSO-3C yield is weaker, other sources provide effective forecast skill (MED 1970(MED -1980;;TNA/Tmax 1990-2005).

For OM, a non-stationary relationship between climate predictors and simulated yield is also seen (Fig. 10b), with much higher correlation values over the 1980s and 1990s (up to 0.67) than in the first period of the 21st century (down to 0.03). These changes can be attributed to the time evolving links between MED/Tmin and OM yield (Fig. 10d).

The results outlined in this sub-section could be associated with multi-decadal changes in the shape and location of teleconnections [START_REF] Jung | Characteristics of the recent eastward shift of interannual NAO variability[END_REF][START_REF] Gómara | Abrupt transitions in the NAO control of explosive North Atlantic cyclone development[END_REF]. For instance, the impact of ENSO on European climate has been shown to be non-stationary in time during the 20th century [START_REF] Mariotti | Euro-Mediterranean rainfall and ENSO-a seasonally varying relationship[END_REF]Greatbach et al., 2004;[START_REF] Rodríguez-Fonseca | A review of ENSO influence on the North Atlantic. A non-stationary signal[END_REF], as is also the case of ENSO impact on Iberian Peninsula crops (Capa-Morocho et al., 2016a,b). An extended assessment considering temporal windows of enhanced prediction opportunity and their underlying factors may require a more complete analysis at this point, which is beyond the scope of this paper. Therefore, it is left as an interesting path for future research.

Final remarks on seasonal forage yield predictions

The results presented here provide support to a future implementation of an in-season forage productivity prediction system at French and European levels, which could be of interest to governments, administrations, chambers of commerce and insurance companies (JPI-Climate, 2019). Grassland production systems are complex and the relationship between forage production and the associated environmental factors are difficult to examine. However, biogeochemical grassland models are ideal for simulating forage production due to their ability to represent the non-linear response of the modelled output to changes of soil, climate and management conditions. To our knowledge, this is the first modeling study to explore the connections between large-scale oceanic-atmospheric circulation and local forage production in Europe. By mobilizing a suite of simulations, it has provided an improved understanding of the response of grassland biomass production to purely climatic patterns. With focus on the French Massif Central, the modeling approach used in this study suggests that: (1) seasonal prediction of grassland productivity from climatic predictors is possible several months prior to cutting; and (2) the links between climate predictors and grassland yield are also not stationary in time.

The relationship between site-specific forage production and NAO and SST patterns identified in this study warrants continue further monitoring to substantiate the role of atmospheric and oceanic changes on grassland dynamics beyond management changes. Further evidence supporting this relationship would permit reasonable inferences about forage production in the form of seasonal (1-5 month) forecasts.

While suggesting the possibility of using grassland models in support of Climate Services, the present study also shows limitations. Our findings are based purely on simulated data to avoid the noise from other factors (e.g., pests, technical development); therefore, they lack evaluation against measured outputs. This is consistent with previous studies that used observed data only for model calibration and validation [START_REF] Capa-Morocho | Crop yield as a bioclimatic index of El Niño impact in Europe: crop forecast implications[END_REF]2016a,b), as simulated production data are considered as the most appropriate to avoid masking plantweather interactions by non-related factors. Moreover, the model used in this study is only a possible realization of alternative approaches to grassland modeling. However, we think that the 'European' calibration/ validation version of PaSim applied for the study site give us, at a minimum, confidence on the reliability of forage production estimates [START_REF] Ma | Regional-scale analysis of carbon and water cycles on managed grassland systems[END_REF]. Another limitation is that grazing schemes were not assessed; this first analysis could only simulate mowing, even though grazing is common in the area. The methodological simplifications adopted in the research (e.g., simulation of potential forage yield, mowing practices) were to some extent necessary owing to the complex nature of grassland systems, as not a single harvest event had to be forecasted, but a continuous biomass production, which in turn, may follow different paths depending on management. Our methodology enables progress towards operative forecasts of forage production. Since this is the first study of this type, it provides ground for further developments, and we expect that a great deal more research will be necessary to unravel the complex relationships at the interface between climate, vegetation and management at different timescales.

Conclusions

The potential links between climate variability, atmospheric and oceanic teleconnection patterns and forage production of a mown permanent grassland system in the French Massif Central (Laqueuille) were analyzed in this study. Following a biogeochemical grassland modeling approach, evidence was provided that skillful seasonal forecasts of forage productivity at the site are achievable using oceanic and atmospheric predictors (e.g., Sea Surface Temperature indices).

This paper provides essential elements to implement a timely and accurate grassland forecasting system in France and, after substantiation in other grassland systems (either with mowing, grazing or both), through gridded simulations over Europe. Such a tool, together with monitoring services (e.g., Copernicus Sentinel satellites; ESA, 2019) and farm economic models [START_REF] Britz | A Highly Detailed Template Model for Dynamic Optimization of Farms -FARMDYN[END_REF], would facilitate to increase on-farm seasonal managerial flexibility, mitigate risks and adopt costeffective local and macro-economic plans (cf. [START_REF] Chen | The signature of sea surface temperature anomalies on the dynamics of semiarid grassland productivity[END_REF]. For instance, by proactively adapting cutting or grazing schedules, stocking rates, irrigation and fertilization plans or production and fallow schemes.

Declaration of Competing Interest

None.

Fig. 1 .

 1 Fig. 1. (a) Location of INRA long-term observational grassland site at Laqueuille (Puy-de-Dôme) and its terrain elevation (in m). Source: NOAA NGDC GLOBE: Gridded 1-km. (b) Mean monthly values of surface solar incident radiation (dark gray -10 W m -2 ), maximum and minimum daily temperatures (red/blue -°C), daily precipitation (cyan -mm day -1 ) and wind speed (magenta -m s -1 ) from the in situ meteorological station (solid lines and bars) and SAFRAN closest grid cell data to Laqueuille (dashed lines and stems). Period 1996-2015 except for solar radiation (2008-2015).

Fig. 2 .

 2 Fig. 2. (a) Regression maps of NAO (shadings: 95% confidence interval; t-test) and mspl-index (contours) on interannual mslp anomalies north of 20°N (hPa SD -1 ). (b) Same as (a) but for Niño 3.4, TNA and MED indices on interannual SST anomalies (°C SD -1 ; red/blue shadings) and AMO index on multi-decadal (> 33-yr; magenta/green shadings). The boxes indicate the SST areas considered for the indices. All shadings indicate 99% conf. int. (t-test).

Fig. 3 .

 3 Fig. 3. (a) Grassland management in MSB simulation (example of year 1982). Black solid line: Potential plant shoot biomass (t DM ha -1 ). Blue solid line: Cumulative fertilization year-round (10 -1 t N ha -1 ). (b) Same as (a) but for 2C management. Green solid line: potential forage yield accumulation (t DM ha -1 ). (c) Same as (b) but for 3C management. (d) Same as (b) but for optimal management.

Fig. 4 .

 4 Fig. 4. (a) Time evolution of annual potential maximum shoot biomass and forage yield (2C, 3C, 3CI and OM) for PaSim simulations (t DM ha -1 ). Period 1959-2015.(b) Time evolution of annual averaged standardized anomalies of surface radiation (srad; W m -2 SD -1 ), Tmax/Tmin (°C SD -1 ), precipitation (Pr; mm day -1 SD -1 ), CO 2 concentration (ppm SD -1 ) and potential maximum shoot biomass (MSB simulation; in t DM ha -1 SD -1 ). (c) Same as (b) but for low-pass filtered (> 11-yr) anomalies of SAFRAN srad, Tmax and NCEP cloud cover (ccover; in % SD -1 ) averaged over mainland France. The multi-decadal NAO and AMO indices (> 33-yr) appear also overlaid.

Fig. 5 .

 5 Fig. 5.(a) Correlation values between high-pass filtered monthly anomalies (12 values per year) of shoot growth rate (gsh; t ha -1 day -1 ) and meteorological variables at the grassland site: srad (W m -2 ), Tmax/Tmin (°C) and Pr (mm day -1 ). 90/95% confidence intervals (t-test) marked with small/big circles. Horizontal lines provide an estimate of the range of non-statistically significant correlations. Vertical lines indicate the starting (dashed) and end (solid) period of maximum shoot biomass measurements (MSB simulation). The mean growth rate values for each month (gsh; t ha -1 day -1 ) are provided in the dotted black line. (b) Same as (a) but for correlations between annual averaged potential yield (MSB; one value per year) and monthly climatic anomalies (12 values per year). (c)-(d) Same as (a)-(b) but for the 3C simulation. Vertical solid lines denote fixed cutting events in time. (e)-(f) Same as (a)-(b) but for the OM simulation. The vertical lines indicate the starting (dashed) and end (solid) period of automatic OM cutting events.

Fig. 6 .

 6 Fig. 6. Correlation values between Laqueuille high-pass filtered potential forage yield (3C simulation) timeseries and January monthly averaged anomalies of: (a) SAFRAN Tmin (shadings and contours; °C). 95% conf. int. in stippling. (b) SAFRAN srad (W m -2 ). (c) NCEP mslp (contours; hPa), HadSST SSTs (shadings over the seas; °C), SAFRAN Tmax (shadings over France) and SAFRAN srad (black contours over France -same as in (b)). 95% conf. int. in all shadings, red/blue thick contours and all black contours. (d)-(l) Same as (a)-(c) but for February, April and September. Period 1959-2015.

Fig. 7 .

 7 Fig. 7. Same as Fig. 6 but for the OM simulation.

Fig. 8 .

 8 Fig. 8. Correlation maps between high-pass filtered 3C potential forage yield (1 data point per year) and January monthly averaged anomalies of SST (HadSST; 95% conf. int. in shadings) and mslp (95% conf. int. in thick contours) from (a) 1959-2015 NCEP, (b) 1959-2001 ERA-40 and (c) 1979-2015 ERA-Interim. (d)-(f) Same as (a)-(c) but for February. (g)-(i) Same as (a)-(c) but for May.

Fig. 9 .

 9 Fig. 9. (a) Same as Fig. 5d but for correlations between TNA, Niño3.4, MED, NAO and mspl-index monthly indices and interannual forage 3C yield anomalies. months: January year -1 to December year 0 (harvest year). (b) Same as (a) but for OM yield anomalies. (c) Same as (a) but for lead/lag correlations of February mspl-index with itself, NAO, MED and Tmin/Tmax/srad/Pr and 3C shoot growth rate (gsh) monthly anomalies at Laqueuille from January to June year 0. (d) Same as Fig 6c but for correlations with February msplindex instead of 3C annual yield. Arrows indicate correlations above 0.1 for positive moisture advection at 850 hPa (NCEP; m g s -1 kg -1 ). (e)-(f) Same as (c)-(d) but for February MED index.

Fig. 10 .

 10 Fig. 10. (a) Stepwise regression model with leave-one-out cross-validation between potential 3C yield anomalies (predictand: bars in t DM ha -1 SD -1 ) and selected predictors (black line) until May of harvest year [TNA September y-1, mspl-index January y0, mspl-index February y0, MED February y0, Tmax Laqueuille April y0 and Niño3.4 May y0]. (b) Same as (a) but for OM. Predictors are selected until July of harvest year [MED February y0 and Tmin Laqueuille July y0]. The 21-yr centered running correlations between predictand/predicted timeseries appear also overlaid (blue lines: 95% confidence interval with filled circles). (c) 21-yr centered running correlations between identified climatic model predictors and 3C yield. Filled circles indicate 95% conf. int. (t-test). (d) Same as (c) but for OM simulation.

Table 2

 2 Summary of PaSim simulations, including information on N fertilization applications, irrigation, number of cuts per year, cutting or maximum shoot biomass dates (day of year), mean potential forage yield and associated standard deviation (in t DM ha -1 ). Weather input in all simulations is SAFRAN for the period 1959-2015.

	Acronym Description
	Simulation

  Dates of cutting events are fixed for all years on days 150 (30 May), 225 (13 August) and 300 (27 October) to try better capture seasonal effects in forage growth (e.g., summer hydrological stress). Potential forage yield value is retained. Fertilization events take place on days 14, 120 and after the two first forage cuts(Fig. 3c -year 1982).

	1.4	2.1	2.2
	8.45	Three annual forage cuts triggered. 8.51	Same as 3C with optimum irrigation conditions. The optimal irrigation method relies on a water stress index from PaSim. 8.68
		3C	3CI
		3 Cuts	3 Cuts + Irrigation	Optimal Management

Acknowledgments

We thank the National Centers for Environmental Prediction, the European Center for Medium-Range Weather Forecasts and Météo-France/Hymex/MISTRALS for providing the NCEP, ERA-40, ERA-Interim and SAFRAN re-analyses. We thank the UK Met-Office Hadley Center for the HadSST database. This study was supported by the project MACSUR -Modeling European Agriculture with Climate Change for food Security (FACCE-JPI), funded by Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA), the metaprogram ACCAF (Adaptation of agriculture and forests to climate change) of the French National Institute for Agricultural Research (INRA) and the project PRE4CAST (CGL2017-86415-R) of the Spanish Ministry of Economy and Competitiveness (MINECO). We also thank the Laqueuille site as part of the SOERE-ACBB project (http://www. soere-acbb.com), funded by the French National Infrastructure (https:// www.anaee-france.fr). Iñigo Gómara was supported by MINECO (Juan de la Cierva-Formación contract; FJCI-2015-23874) and Universidad Politécnica de Madrid (Programa Propio -Retención de Talento Doctor). Iñigo Gómara's research stay at INRA UREP (September-December 2017) was also funded by UPM (Programa Propio -Ayudas al personal docente e investigador para estancias breves en el extranjero). We would like to thank Olivier Darsonville, Rémi Perrone, Iris Lochon, Katja Klumpp, Catherine Picon-Cochard (INRA UREP), Pere Quintana-Seguí (Observatori de l'Ebre) and Roberto Suárez-Moreno (Columbia University) for their helpful cooperation during the progress of this study. Finally, we are indebted to the two anonymous reviewers, whose pertinent comments and suggestions have contributed to substantially improve this manuscript.

Supplementary material

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.agrformet.2019.107768.