Climate Change Effects upon Pasture in the Alps: The Case of Valtellina Valley, Italy

: In this study, we assessed the potential effects of climate change upon the productivity of mountain pastures in the Valtellina valley of Italy. Two species, Trisetum ﬂavescens and Nardus stricta , among the most abundant in Italian pastures, were chosen for the simulation of low-and high-altitude pastures, respectively. We introduced some agroclimatic indices, related to growing season parameters, climate, and water availability, to evaluate the impacts of climate change upon pasture production. First, the dynamic of the pasture species was evaluated for the present period using the climate-driven, hydrologically based model Poli-Hydro , nesting the Poli-Pasture module simulating plants growth. Poli-Pasture was validated against yield data, at province scale, and at local scale. Then, agroclimatic indices were calculated. Subsequently, IPCC scenarios of the Fifth and Sixth Assessment Reports (AR5 and AR6) were used to project species production and agroclimatic indices until the end of the 21st century. In response to increased temperature under all scenarios, a large potential for an increased growing season length and species yield overall (between +30% and +180% for AR5 at 2100) was found. Potential for decreased yield (until − 31% for AR5) is seen below 1100 m asl in response to heat stress; however, it is compensated by a large increase higher up (between +50% and +140% for AR5 above 2000 m asl). Larger evapotranspiration is foreseen and larger water demand expected. However, speciﬁc (for hectares of pasture) water use would decrease visibly, and no signiﬁcant water limitations would be seen. Results provide preliminary evidence of potential livestock, and thereby economic development in the valley at higher altitudes than now.


Introduction
In mountain areas, pastures and farming systems are paramount important activities for local communities, a source of income for local development, and a key feature of local ecosystems dynamics [1]. Pasture management has positive effects on land sustainability, maintaining the landscape and cultural value and supporting biodiversity and soil fertility, thereby reducing soil loss and natural risks. On the contrary, land abandonment brings the growth of shrubs, grassland biodiversity loss, and increased erosion, wildfire, and avalanche risk, as well as the loss of traditional productive activities and typical landscapes [2]. During 1990-2010, ca. 17% of the Italian alpine pastures were abandoned [2].
Projected trends until 2100 under the IPCC scenarios depict a large decrease in snowfall and snow cover during winter [9] (high confidence [3]), with ever earlier melting, larger The study area covers a large altitude range, from 225 to 3610 m asl. Valtellina has a temperate climate, with no dry seasons and a warm summer, classified as Cfb in the Köppen-Geiger climate classification [43]. Yearly precipitation P in the area averages 1400 mm y −1 (1422 mm y −1 below 2000 m asl, 1383 mm y −1 above 2000 m asl). Precipitation has no large lapse rate in elevation. Usually, however, P is more abundant in the western part of the valley (nearby Como Lake, e.g., Morbegno, 1137 mm y −1 ) than in the eastern part (e.g., Livigno-Foscagno, 670 mm y −1 , see also [44]). The mean annual temperature, T, is +1.5 °C above 2000 m asl and +5.0 °C below 2000 m asl. It reaches 0 °C at Livigno La Vallaccia (2660 m asl) and +14 °C in Morbegno (230 m asl). The summer temperature peaks in July at around +23 °C in Morbegno and +8 °C in Livigno. The winter temperature reaches −8 °C and +3 °C, respectively.
Pasture and grassland cover an area of 522 km 2 within the basin (Figure 1). A large spatial variability of pasture productivity is present because of differences in temperature (i.e., altitude), soil fertility, precipitation, and soil moisture conditions [45]. Usually, productivity (tons of dry product per hectare per year) may range within 0.5-6.5 td ha −1 y −1 . In Valtellina, bovines, ovines, and caprines graze on grasslands and pastures. In particular, during the spring and early summer, animals graze at low altitudes, and then during the summer grazing is moved upward to the pasture lands. The grazing may be free or controlled. In the first case, animals, especially bovines, can move for preferred grass. The second case applies principally to ovines/caprines and animals gathered in defined The study area covers a large altitude range, from 225 to 3610 m asl. Valtellina has a temperate climate, with no dry seasons and a warm summer, classified as Cfb in the Köppen-Geiger climate classification [43]. Yearly precipitation P in the area averages 1400 mm y −1 (1422 mm y −1 below 2000 m asl, 1383 mm y −1 above 2000 m asl). Precipitation has no large lapse rate in elevation. Usually, however, P is more abundant in the western part of the valley (nearby Como Lake, e.g., Morbegno, 1137 mm y −1 ) than in the eastern part (e.g., Livigno-Foscagno, 670 mm y −1 , see also [44]). The mean annual temperature, T, is +1.5 • C above 2000 m asl and +5.0 • C below 2000 m asl. It reaches 0 • C at Livigno La Vallaccia (2660 m asl) and +14 • C in Morbegno (230 m asl). The summer temperature peaks in July at around +23 • C in Morbegno and +8 • C in Livigno. The winter temperature reaches −8 • C and +3 • C, respectively.
Pasture and grassland cover an area of 522 km 2 within the basin (Figure 1). A large spatial variability of pasture productivity is present because of differences in temperature (i.e., altitude), soil fertility, precipitation, and soil moisture conditions [45]. Usually, productivity (tons of dry product per hectare per year) may range within 0.5-6.5 t d ha −1 y −1 . In Valtellina, bovines, ovines, and caprines graze on grasslands and pastures. In particular, during the spring and early summer, animals graze at low altitudes, and then during the summer grazing is moved upward to the pasture lands. The grazing may be free or controlled. In the first case, animals, especially bovines, can move for preferred grass. The second case applies principally to ovines/caprines and animals gathered in defined areas, where grass is at the optimal growth stage. Animals optimize their consumption, and it is possible to control the growth of good pasture species [46,47].

Data and Methods
In Figure 2, two flow charts of the methodology used here are reported, while in Appendix A, a summary of the two models used here is given. simulated availability of water for plant growth.
The principal purposes of the work are: (i) to simulate distributed pasture productivity in a watershed to have results on a large area and to identify differences in zones; (ii) to study the timing of plant growth, using both a simulation with both fixed growing season and variable growing season; (iii) to understand if synthetic indices could be useful to identify causes of growth limitation.
After the calibration of the model with the two growing season schemes, the agroclimatic indices were defined. To conclude, future projections were performed using scenarios of the IPCC AR5 and AR6. It was possible to analyze different results of the two assessment reports, the overall impact of climate change on pasture dynamic and growth, and the efficacy of the chosen agroclimatic indices. Poli-Hydro (hydrology) and Poli-Pasture (pasture growth) components and interconnections are reported in the frame on the top, named "present". Generation of future projections is presented in the frame below, called "future". T, P, and Rad are meteorological data: temperature, precipitation, and radiation, respectively. S(t) is the soil water content, and ETeff is actual evapotranspiration. (b) Detailed flow chart of Poli-Pasture model: depicting operation at a generic time step t. Inputs are data, Poli-Hydro-simulated variables at time step t, and Poli-Pasture-simulated variables at the previous time step. Outputs are used for the agroclimatic indices calculation and as input to Poli-Hydro and Poli-Pasture at the following time step.

Hydrological and Pasture Model
The Poli-Hydro model is essentially a spatially semi-distributed (i.e., with distributed budget, simplified cell-by-cell flow routing, independent cells, and lateral flows neglected) model, which applies mass conservation (continuity) equations to water content in the soil between two consequent time steps. The model (see [48,49]) is suitable for assessment of water budget and cryospheric dynamics in high-altitude areas (e.g., [50]). The study area was partitioned in a grid with cells of 1 × 1 km 2 , and for each time step (here, 1 Poli-Hydro (hydrology) and Poli-Pasture (pasture growth) components and interconnections are reported in the frame on the top, named "present". Generation of future projections is presented in the frame below, called "future". T, P, and Rad are meteorological data: temperature, precipitation, and radiation, respectively. S(t) is the soil water content, and ETeff is actual evapotranspiration. (b) Detailed flow chart of Poli-Pasture model: depicting operation at a generic time step t. Inputs are data, Poli-Hydrosimulated variables at time step t, and Poli-Pasture-simulated variables at the previous time step. Outputs are used for the agroclimatic indices calculation and as input to Poli-Hydro and Poli-Pasture at the following time step.
In practice, the hydrological model Poli-Hydro is coupled with a pasture module Poli-Pasture to be able to simulate pasture dynamics starting from some outputs of the hydrological model: e.g., the spatialized daily temperature, the on-ground radiation, and the simulated availability of water for plant growth.
The principal purposes of the work are: (i) to simulate distributed pasture productivity in a watershed to have results on a large area and to identify differences in zones; (ii) to study the timing of plant growth, using both a simulation with both fixed growing season and variable growing season; (iii) to understand if synthetic indices could be useful to identify causes of growth limitation.
After the calibration of the model with the two growing season schemes, the agroclimatic indices were defined. To conclude, future projections were performed using scenarios of the IPCC AR5 and AR6. It was possible to analyze different results of the two assessment reports, the overall impact of climate change on pasture dynamic and growth, and the efficacy of the chosen agroclimatic indices.

Hydrological and Pasture Model
The Poli-Hydro model is essentially a spatially semi-distributed (i.e., with distributed budget, simplified cell-by-cell flow routing, independent cells, and lateral flows neglected) model, which applies mass conservation (continuity) equations to water content in the soil between two consequent time steps. The model (see [48,49]) is suitable for assessment of water budget and cryospheric dynamics in high-altitude areas (e.g., [50]). The study area was partitioned in a grid with cells of 1 × 1 km 2 , and for each time step (here, 1 day) the required variables were evaluated. The Poli-Pasture module [24,33,51] iteratively receives from Poli-Hydro the daily value of soil water content, to evaluate the growth of vegetal/agricultural species for the same time steps, and cells. The Poli-Pasture module returns to Poli-Hydro a daily value of leaf area index (LAI), used to assess correct daily evapotranspiration, fraction of soil with vegetal coverage, and soil water content (used for the vegetation growth).
Potential evapotranspiration ET max was calculated using Priestley-Taylor's formula in pasture areas and Hargreaves' formula in non-pasture areas. Then, the model estimated actual (soil moisture-driven) evapotranspiration ET eff . Comparing the two values (ET max and ET eff ), it was possible to highlight water stress resulting from low precipitation and/or high temperature. Evapotranspiration, in the presence of vegetation, depends upon the LAI, which was calculated daily, based upon the daily vegetation growth, its vegetative stage, and soil water content [52].
The growth stages of the pasture species were assessed based upon accumulation of thermal time stages (degree-days DD [ • C]). During the growing season (GS) [30], a new stage is reached when the accumulated thermal time reaches a proper threshold. If the daily temperature is below a base temperature T base , no thermal time is accumulated, and if the daily temperature is larger than a cutoff temperature T cutoff , the growth is limited. Model tuning was attained by varying degree-day thresholds for each growth stage (beginning of growth, flowering, maturity, end of growth). The Poli-Pasture module estimates daily production of biomass for each species, i.e., the minimum value between (i) water-dependent growth G TR and (ii) solar radiation growth G R: For the simulation, proper dates for the growing season start/end needed to be chosen. GS started on April 1st in LowAlt and on April 30th in HighAlt, while the end of GS was taken on September 30th in both zones [53]. At low altitude, the Trisetum flavescens is harvested just before the maturity, while Nardus stricta at high altitudes is fed to animals just after the flowering, as reported in literature (e.g., [45,54]). After cut/feeding, growth starts again, and the cycle is repeated until the end of the GS. The total annual production in a cell is calculated as the sum of peak biomass in each growth cycle.
Then, the same simulation was carried out considering a variable growing season. GS would start when the mean temperature would be higher than T base (for each species) for 10 consecutive days, and it would stop when T would be lower than the mean temperature of September 30th (during the period 2006-2019) (i.e., the expected end of the GS) for 10 consecutive days. This may account for yearly potential variations of the growing season length, in response to specific climate conditions. In doing so, it was possible to (i) improve the results of simulation, and (ii) evaluate the impact of the future modified climate upon species' growing season length, therefore and productivity.

Weather Data
The Poli-Hydro model is driven by meteorological data, namely temperature, (total) precipitation (then properly split into solid/liquid), and radiation, and it gives as output estimates the soil moisture, evapotranspiration, and (total, overland plus sub-superficial) runoff. It was run here at the daily scale, suitable for pasture simulation (e.g., [24]).
Weather data are made available by ARPA Lombardia (the regional agency for environmental protection, www.arpalombardia.it, accessed on 3 November 2022). For this study, 15 ARPA stations were used (see Figure 1 and Table 1). During 2003-2019, each station measured daily data of (air) temperature T. A total of 8 stations provide total precipitation P (heated rain gages, including solid precipitation), while 7 others separately provide rainfall R and snow depth H (the latter via sonic snow gauging), always at the daily scale. When snow depth H data are available, one can assess daily new snow depth H n and then convert it into new snow water equivalent (SWE n ). Here we took a mean snow density of 115 kg/m 3 , valid in the area [55,56]. Table 1. Meteorological stations considered for the model: elevation, position, and measured variables. T is temperature, P is total precipitation (heated rain gauge), R is (liquid) rainfall, H is snow depth. The location of the stations is visible in Figure 1.  The Poli-Hydro model needs distributed daily meteorological data, i.e., a value of temperature and precipitation for each time step and each cell. Thus, AWS daily data were spatialized considering Thiessen polygons around each AWS and applying a monthly altitudinal gradient (for T and P) calculated from the AWS data and their elevation. Solar radiation was calculated on a distributed basis (i.e., in every cell of the model), based upon theoretical extraterrestrial radiation, topographically corrected for shading ( [57]).

Land Cover and Soil Properties
A land cover map was used to estimate the maximum soil retention potential S Max (necessary to pursue Poli-Hydro simulation) within the study area, based upon the SCS-CN method [58], and also to assess the area covered with vegetation (pastures in particular) and with rocks and ice, etc. This map was made available under the umbrella of the Corine Land Cover (CLC, 2012) experiment of the European Copernicus Programme. The Poli-Hydro model requires as inputs the hydraulic properties of soil, namely wilting point θ w , field capacity θ l , saturation θ s , hydraulic conductivity K, in turn depending on soil texture [59] and the depth of the active soil layer, considered here constant for the entire area, based on the average properties of the Valtellina territory. Soil properties used for hydrological modeling were also made available in fulfilment of other studies of the area of Valtellina [48,57,60].

Pasture Productivity
The Poli-Pasture module runs jointly with the Poli-Hydro model to calculate pasture species productivity. Poli-Pasture was tuned here using pasture yield statistics, as reported by ISTAT (i.e., the National Institute of Statistics), aggregated for each province of Italy. Valtellina and Valchiavenna (included in the study area) cover the whole territory of Sondrio province, and accordingly the total production simulated by the model in the considered area is comparable to the data given by ISTAT for that province.
Pasture lands in Valtellina are substantially divided into two areas, namely highaltitude pastures, above 2000 m asl (henceforth pasture zone HighAlt), and grasslands at low altitude, below 2000 m asl (pasture zone LowAlt).
In fulfilment of the MOUPA project, data of pasture (Nardus stricta) biomass and relative species abundance were collected in the Dosdè area ( Figure 1) during summer 2019 (2 different sites, Dosdè 1, 2100 m asl, and Dosdè 2, 2500 m asl, Confalonieri R., personal communication, 2020). Based upon the available literature, and upon results from collected samples under the umbrella of the MOUPA project, within the HighAlt/LowAlt zones we could identify some more abundant species. Above 2000 m asl, the most abundant species retrieved was Nardus stricta (henceforth Nar.), while below 2000 m asl Trisetum flavescens (henceforth Tris.) was predominant [53,54]. In particular, in the first considered site of Dosdè, Nardus stricta is the second most abundant with a presence of 10.25% in July and 34.47% at the end of August; in the second site, it was the most abundant species with a percentage of 84.45% in July and 77.44% in August. Therefore, these species were preliminarily considered as representative of the pasture vegetation in the area [45,61] and their growth was simulated. More realistic pasture growth simulation may require modeling of multispecies growth and their intraspecies and interspecies competition [25,62]. However, for our purpose of preliminarily assessing potential impacts of climate upon pasture productivity in this area, the hypothesis of monoculture pastures may be seen as representative of bulk pasture behavior (e.g., [24]). For the simulation, specific parameters of growth were collected for both species [52,[63][64][65]. In Table 2, the parameters used for simulation with the Poli-Pasture model are reported for both species. The values of these parameters were mostly collected from literature (see e.g., [24,26,33,[66][67][68]), and some were tuned with manual calibration. The reported values of DD factors were cumulated since the beginning of the GS. Table 2. Poli-Pasture module parameters: their value for Nardus stricta and Trisetum flavescens. Parameters in bold calibrated against yield data from ISTAT. Other parameters taken based upon literature. a: [68], b: [67], c: [33], d: [24], e: [26]. DD factors are taken since growth season GS start. ISTAT reported estimates of total annual production Y y [t] (2006-2017) in the Sondrio province, but pasture area therein varied during the years. Indeed, the area used for grasslands and pastures decreased during the study period, according to ISTAT information. Such changes are not justified in the ISTAT report and might have been influenced by social and economic factors, not amenable to modeling here, such as land abandonment for personal reasons, generational changes, or lack of public funding. Accordingly, specific (to area) production Y was then used for calibration purposes. Model tuning was pursued according to an objective score, namely the percentage mean error Bias. An acceptable calibration would be attained if |Bias|<10%. Given the small amount of data for tuning (i.e., 14 yearly values of Y), other more sophisticated adaptation statistics (root-mean-square error RMSE and root-mean-squared percentage error RMSE % ) were not used for tuning, but were, however, calculated and are reported further on.

Variable
It was not possible here to verify growth timing and spatial distribution of biomass in the area, given the lack of in situ data or spatially distributed information. However, the growth patterns of the two species were compared against those found in other studies about pasture [24,45], and Y was compared against literature values. At high altitudes in the Alps, productivity of pastures was estimated to be near to 3-4 t ha −1 , while below 1000 m asl higher values could be found, around 8 t ha −1 [45,54]. Here, it was also possible to gather an in situ estimation of specific annual production (of Nardus stricta) at high altitudes in the area of Valtellina valley. This came from the available literature, specifically in Alpe Boron (Figure 1), during 2003 and 2004 [69], and from production data collected in the Dosdè area in fulfilment of the MOUPA project, and could be used here for the benchmark.

Climate Projections and Future Simulations
Once properly tuned, Poli-Hydro and Poli-Pasture were used to project pasture production during 2020-2100 (i.e., after the end of the present period 2003-2019). After a preliminary screening, 21 scenarios were generated, combining different GCMs and RCPs/SSPs, to analyze the response of our pasture species growth to different pathways of the climatic drivers. Particularly, the selected GCMs are able to acceptably describe the climate of Northern Italy, especially in terms of the seasonality of precipitation [18,29,70]. However, local downscaling is required to properly mimic local weather statistics at specific points [17].
The simulation was conducted under a variable GS, under three representative concentration pathways (RCPs) of the Fifth Assessment Report of IPCC, RCP 2.6, RCP 4.5, RCP 8.5. Three atmospheric and ocean general circulation models (AOGCMs) were used, downscaled to the area of study, from the Coupled Model Intercomparison Project, release 5 [ [78]).
In general, GCMs display a variable spatial resolution, in the order of 100 × 100 km 2 or so, and downscaling needs to be pursued to gather daily series of temperature and precipitation at the spatial resolution of 1 × 1 km 2 for the Poli-Hydro model. For precipitation, a stochastic space random cascade model (SSRC, [70]) was used. This method provides downscaled series that are statistically consistent (i.e., displaying the same properties of mean, variance, and intermittence for rainfall) with the locally observed data, but conserve the trends as provided by the GCM models. Temperature series were downscaled using a ∆T approach [17], in which an average monthly difference of temperature between the GCM series and the observed series calculated in the control run CR period (2006-2017) was applied to the projected GCM series, largely depending upon local altitude. By downscaling each GCM and scenario separately, spatially distributed, daily projected precipitation and temperature scenarios were obtained and used to feed the Poli-Hydro + Poli-Pasture model. Model parameters are held constant in the future in lack of any other hypothesis, and the area of pasture is the same, although chances for changes may occur.
It was decided to simulate both AR5 and AR6 because of the differences in their projections. For example, AR6 shows more pronounced increases of temperature with respect to AR5 for all GCMs and all scenarios. Considering precipitation, GCMs for the AR6 seem to be more unanimous in the projections [79]. Generally, GCMs for the AR5 project a large increase or large decrease in precipitation, while for the AR6 all project a contained reduction.
Heat waves frequency (number of consecutive days with T > T cutoff ) Number of days in GS with precipitation > 10 mm d 10 d f, g AI4 Total precipitation in growing season P GS mm g AI5 Annual species productivity Y t a, c, e AI6 ET efficiency (in GS) ET eff /ET max,GS mm/mm -AI7 ET relative (in GS) ET eff /P GS mm/mm -AI8 Specific (green) water footprint (in GS) Index AI1 provides the length of the growing season, and clearly, it was calculated only under the variable GS mode. Index AI2, or frequency of heat waves, and AI3, days with heavy precipitation during the growing season GS, represent the variability of stress factors for the species. Index AI4 of total precipitation indicates abundance/lack of rainfall for vegetation growth. Index AI5 provides an indication of (total) species biomass in the area. AI6, namely the ratio of yearly actual-to-potential evapotranspiration in the GS season, assesses the efficiency of species water use, in turn depending upon water availability and distribution in time (e.g., [34]). If this ratio is close to 1, plants use most of the available water, with large evapotranspiration efficiency. In AI7, relative ET is an indication of the necessary evapotranspiration with respect to the available water (precipitation, i.e., rainfall) in the GS period. This depends upon temperature (driving ET max,GS ) and available precipitation P. Specific water footprint AI8 is an indication of how much water (ET eff ) is needed to be used for the production of a ton of biomass (e.g., [80]).
AI values were calculated under present and future conditions for each RCP/GCM under AR5/AR6. The indices AI1-AI4 could be assessed using only climate data. Accordingly, climate data from the available stations were used for calculation under present conditions, while climate projections were used for future conditions. The indices AI5-AI8 explicitly depict the pasture species performance, and they were calculated using Poli-Pasture model's outputs under present and projected climate. Poli-Pasture calibration was pursued by iteratively changing the tuning parameters (Table 2) to mimic the observed values of specific annual production Y from ISTAT (2006-2017). Clearly, the robustness of the parameters' tuning is of interest, and a sensitivity analysis within the parameters' viable space could be pursued [81]. However, here we were more interested in choosing the parameters' values that were most suitable to properly simulate the plants' dynamics. This was reasonable, given also that we possessed some data, useful to reasonably constrain the uncertainty in model tuning. A preliminary analysis here, and in other former studies, demonstrated that DD factors are in practice most affecting the performance of the model; thus, proper tuning is required.
The calculated percentage error, or Bias % , for the simulation with fixed GS was +10.82%, considering the production per hectare Y (−4.90% considering the total production Y y ). In the simulation with variable GS, Bias % was −0.18% for Y (−9.93% for Y y ). Such a difference between Y, and Y y may depend on changes in the pasture areas during the study period, as reported above. Goodness-of-fit statistics for Y for Poli-Pasture-i.e., Bias % , random mean square error, absolute RMSE, and percentage RMSE % -are reported in Table 4. Determination coefficient R 2 was not calculated. Due to (possibly unlikely) constant values of Y, as declared by ISTAT during 2006-2007 and 2008-2012, the variance of the observed sample is very low, and accordingly, calculation of R 2 makes little sense. In Figure 3, observed and simulated values of Y in the CR period are reported. Poli-Pasture calibration was pursued by iteratively changing the tuning parameters ( Table 2) to mimic the observed values of specific annual production Y from ISTAT (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Clearly, the robustness of the parameters' tuning is of interest, and a sensitivity analysis within the parameters' viable space could be pursued [82]. However, here we were more interested in choosing the parameters' values that were most suitable to properly simulate the plants' dynamics. This was reasonable, given also that we possessed some data, useful to reasonably constrain the uncertainty in model tuning. A preliminary analysis here, and in other former studies, demonstrated that DD factors are in practice most affecting the performance of the model; thus, proper tuning is required.
The calculated percentage error, or Bias%, for the simulation with fixed GS was +10.82%, considering the production per hectare Y (−4.90% considering the total production Yy). In the simulation with variable GS, Bias% was -0.18% for Y (−9.93% for Yy). Such a difference between Y, and Yy may depend on changes in the pasture areas during the study period, as reported above. Goodness-of-fit statistics for Y for Poli-Pasture-i.e., Bias%, random mean square error, absolute RMSE, and percentage RMSE%-are reported in Table 4. Determination coefficient R 2 was not calculated. Due to (possibly unlikely) constant values of Y, as declared by ISTAT during 2006-2007 and 2008-2012, the variance of the observed sample is very low, and accordingly, calculation of R 2 makes little sense.
In Figure 3, observed and simulated values of Y in the CR period are reported.  For further model validation, in Figure 3 we reported in situ estimates of biomass of Nardus stricta in Alpe Boron (2003)(2004) and Alpe Dosdè (2019, 2 sites), and the estimates from Poli-Pasture model, in the corresponding cells of our 1 × 1 km 2 grid. Visibly, except for the year 2003 when some difference is spotted, in Boron the model simulates a specific productivity of 3.51 t/ha (fixed GS) and 2.92 t/ha (variable GS) instead of the observed value 1.03 t/ha. The Poli-Hydro model matches reasonably well with local (i.e., within a specific cell) estimates of pasture yield. In Table 5, the values of the difference in terms of specific productivity Y are reported for both the sites of Boron (for year 2004) and Dosdè (as average of the two sites). Such data were not used for model tuning, and accordingly, these results indicate an acceptable validation of the model.

Uncertainties in Calibration Data
ISTAT data of pasture productivity give the yearly sum of the total production of three classes of pasture/grassland, namely permanent grasslands and pastures, poor pastures, and other pastures. These data are collected with an estimative methodology, based upon the evaluations of local experts, producers associations, questionnaires, and auxiliary data, etc. The estimates for 2006 and 2007 derive from a different campaign, with different guidelines for the statistics of plant-related products. Moreover, the class "poor pastures" appears only since 2013. Poor pastures have a lower specific production Y with respect to other classes of pasture, and they contribute to the reduction of specific production after 2013.
Because of these uncertainties related to ISTAT data, in addition to the comparison of specific productivity, we also compared the growth rate between the simulation results and the observation data. As reported in Figure 4, looking at the ISTAT data of productivity, only for permanent grassland and other pastures in the period 2006-2017 was there a growth rate of +0.11 t/ha/y, coherent with the result of our simulation (growth rate of +0.07 t/ha/y in the simulation with fixed GS, and of +0.10 t/ha/y in the simulation with variable GS). Moreover, as we can see in Figure 3, in 2015-2020 the overall (considering all three ISTAT classes of pastures) specific production also has a growing trend equal to +0.12 t/ha/y.  The yield of pasture species Y and Yy could also be evaluated as per altitude belts, ranging from 200 to 2600 m asl. The study area was divided into eight belts of 300 m of altitude, each one displaying a different area. In Figure 5, the distribution of Yy, and Y as per altitude belts is reported for the year 2017; however, the altitude pattern is very similar Given the complexity of the vegetation modeling exercise and the potential for some uncertainty in ISTAT estimates of both yield and of pasture-covered areas, Poli-Pasture estimates seem acceptable, at least on average. Furthermore, given the fact that we are interested in assessing relative changes in the yield potential of our pasture species under climate change, the model's estimates can be taken here as useful for the purpose.
The yield of pasture species Y and Y y could also be evaluated as per altitude belts, ranging from 200 to 2600 m asl. The study area was divided into eight belts of 300 m of altitude, each one displaying a different area. In Figure 5, the distribution of Y y , and Y as per altitude belts is reported for the year 2017; however, the altitude pattern is very similar for every year during the present period.
Similar results were found when considering the two main altitude zones but with variations in HighAlt larger than in LowAlt. The range of variation ΔY (mean, min, max) among all models and RCPs as per altitude belts during 2041-2050 and 2091-2100 is reported ( Figure 6). In HighAlt (>2000 m asl), the percentage variation ΔY% is larger than in LowAlt. However, considering the absolute values of Y, future productivity would still be smaller in HighAlt (Nardus stricta) than in LowAlt (Trisetum flavescens). The main difference between AR5 and AR6 is that, according to the latter, during 2041-2050 and 2091-2100 similar changes ΔY would be seen (Figure 6b), differently from AR5 (Figure 6a). Above 1700 m asl or so, AR6 projects a larger variation for 2041-2050 and less for 2091-2100 than AR5. Moreover, in some belts (200-500 m asl and 800-1100 m asl) under AR5, productivity would slightly reduce until -18% on average (Figure 6a, 3rd belt, 2041-2050) within a range reaching −31%, and the reduction would be larger in 2041-2050 than in 2091-2100. In AR6, there would be a reduction for the same belts, slightly larger in 2091-2100 than in 2041-2050, and smaller than in AR5, namely reaching until -27.40% ( Figure  6b, 3rd belt, 2091-2100). On average, however, reduction would only occur in the first belt (−10%).  A largest yield Y can be seen between 1700 m asl and 2000 m asl (with Tris.) and between 2000 and 2300 for Nar., where, however, the specific production Y is quite low (but clearly a large area for pasture cropping is available).

Future Pasture Species Productivity
Given the outputs of the model tuning exercise, we decided to simulate the future growth of our pasture species under a variable GS mode, especially to be able to highlight changes in the duration of the suitable season for species growth. Generally, an increase in productivity is expected. This is principally due to an increase in T. Indeed, temperature increases potential evapotranspiration ET max . If enough water is available, actual evapotranspiration ET eff also increases, thus increasing ET-dependent biomass growth (Equation (A8)). Increased temperatures lead to a quicker achievement of saturation of the degree-day, with consequently faster maturity and possibly more growth cycles during GS. Looking at AR5, for all models, under RCP 2.6 the differences between the period 2041-2050 and the period 2091-2100 are small, but for RCP 8.5, ∆Y is significantly larger for 2091-2100 than for 2041-2050. In general, variations during 2091-2100 are larger for AR6 than for AR5. Similar results were found when considering the two main altitude zones but with variations in HighAlt larger than in LowAlt. The range of variation ∆Y (mean, min, max) among all models and RCPs as per altitude belts during 2041-2050 and 2091-2100 is reported ( Figure 6). In HighAlt (>2000 m asl), the percentage variation ∆Y % is larger than in LowAlt. However, considering the absolute values of Y, future productivity would still be smaller in HighAlt (Nardus stricta) than in LowAlt (Trisetum flavescens). The main difference between AR5 and AR6 is that, according to the latter, during 2041-2050 and 2091-2100 similar changes ∆Y would be seen (Figure 6b), differently from AR5 (Figure 6a). Above 1700 m asl or so, AR6 projects a larger variation for 2041-2050 and less for 2091-2100 than AR5. Moreover, in some belts (200-500 m asl and 800-1100 m asl) under AR5, productivity would slightly reduce until −18% on average (Figure 6a, 3rd belt, 2041-2050) within a range reaching −31%, and the reduction would be larger in 2041-2050 than in 2091-2100. In AR6, there would be a reduction for the same belts, slightly larger in 2091-2100 than in 2041-2050, and smaller than in AR5, namely reaching until −27.40% (Figure 6b, 3rd belt, 2091-2100). On average, however, reduction would only occur in the first belt (−10%).

Agroclimatic Indices
The values of AIs were assessed in every year of the CR period for Trisetum flavescens (<2000 m asl, LowAlt) and Nardus stricta (>2000 m asl, HighAlt). Then, the AIs were calculated for the future projections during 2041-2050 and 2091-2100, for each model and RCP, with variable GS mode. All AIs are reported in Figure 7. The figure presents the average values (among all GCMs) of the indices for each scenario, also benchmarked against values in the control run period CR (2006-2017) for Nardus stricta (HighAlt, >2000 m asl) and Trisetum flavescens (LowAlt, <2000 m asl), at the middle (2041-2050) and at the end (2091-2100) of the 21st century.
Generally, a potential extension of GS (i.e., a longer growing season, AI1, Figure 7a) during 2041-2050, and 2091-2100 is seen with few exceptions under AR5. Moreover, AR6 projections show a larger extension of GS with respect to AR5 projections. Longer GSL is due both to an earlier onset of GS and to a later end of GS due to higher temperatures in the spring and fall. This clearly contributes to the increase in yield Y (AI5, Figure 7e) because of (i) the longer time for growth and (ii) the chances for more growth cycles (harvest/feed). However, a significant increase in temperature, with more frequent heat waves (AI2, Figure 7b), may be a limitation for plants' growth and may justify the potential decrease in Y in Figure 6 (low altitude belts-LowAlt). Comparing projections of AI3 (days with P > 10 mm, Figure 7c) and AI4 (total precipitation during GS, Figure 7d), mostly a reduction of precipitation brings a reduction of intense precipitation days, and vice versa, with some exceptions, both during 2041-2050 and 2091-2100.
Index AI6, ET efficiency, or the ratio ET eff /ET max (Figure 7f) increases in 2041-2050 and 2091-2100 in both HighAlt and LowAlt, in spite of mostly decreasing precipitation. Indeed, in our simulations, ET eff increases more than ET max , possibly in response to (i) the longer growing season and (ii) the decreased incidence of extreme events (AI1, AI4) and subsequently more regular distribution of precipitation/soil moisture.
AI7, relative ET, or ET eff /P GS increases mostly (Figure 7g) in both HighAlt and LowAlt in 2041-2050 and 2091-2100. Under increased temperature, especially in the highest areas, more water from precipitation is used to fulfil ET requirements for plants' growth (resulting in increased yield, as reported in Figure 6).
Specific water footprint AI8, namely ET eff /Y (Figure 7h), may increase or decrease depending on the combination of ET eff , and Y dynamics.
Comparing AR5 and AR6 projections, relatively to AI6 no differences could be seen. However, a larger increase in ET efficiency is expected at high altitude with respect to low altitude (under 2000 m asl).
On the contrary, considering AI7 (relative efficiency), a larger increase is projected by AR6 scenarios, both above and under 2000 m asl (HighAlt and LowAlt). Considering there are no marked differences in the precipitation quantity during GS based on the selected AR, this means a larger increase in ET eff for the AR6.
Considering the water footprint AI8, under the AR5 our GCMs usually project an increase for low altitude (LowAlt) and a decrease for high altitude (HighAlt) for most RCPs (with the exception of the ECHAM6.0 model). On the contrary for AR6 scenarios, the GCMs show a decrease in specific water footprint for both low and high altitude, generally with few exceptions.
In spite of the increase in ET eff , highlighted by AI6 and AI7 in particular for the AR6, and in general for high-altitude areas (both AR5 and AR6), the large increase in productivity at high altitude and in particular for the AR6 explains the reduction of the specific WF (AI8). On the contrary, the decrease in productivity at low altitude explains the reduction of AI8 in LowAlt.

The Calibration of Poli-Hydro+Poli-Pasture Model
The use of the Poli-Pasture model, associated to the hydrological model Poli-Hydro, allowed here to link growth of our pasture species to the climate and water budget of the study area. The Poli-Hydro model is suitable for the simulation of high-altitude watershed hydrology, as in the case of Valtellina and Valchiavenna [29,48,57,60], thanks to the accurate formulation of snow melt and, in general, of water balance components. In [29], the authors report proper calibration of the hydrological model Poli-Hydro for the considered area. Poli-Pasture may therefore be used to simulate growth of high-altitude pasture species, differently from other models used for low grasslands [82][83][84][85][86][87]. Particularly, snow cover dynamics affect the growing season of pasture species and water availability, and simulation is thereby important in our target areas.
With proper hydrological simulation as a background, the modeled pasture species growth seems acceptably representative. The main issue with model tuning here was data availability, given that we could use only bulk pasture yield estimates at the large scale of the Sondrio province. These data carried information from areas of permanent/rich pasture and from other areas of transient/poor pasture with lower yield; however, they were not documented properly. The model performance was influenced by the lack (or paucity) of local data about pasture yield and management of pasture areas, and the large differences in altitude and climatic conditions within the study area (see e.g., [51]).
The results here provide errors smaller than 10% on average. Statistics of random error, i.e., RMSE % , near to 30% seem acceptable [33,88]. However, validation against some local, sparsely available Nardus stricta biomass estimates at high altitudes provides encouraging results. In addition, our findings are coherent with recent studies on pasture biomass ( [24], focusing upon pastures in mountain areas of Sardinia; [26,69], focusing on Valtellina).

Projections of Potential Changes in Pasture Dynamic and Productivity
The results of our future projections seem well aligned with recent literature. GS would be longer (Figure 7a), and production of species would increase (Figure 4). Water-use efficiency would increase, as shown in Figure 7f. An upward shift is seen here in terms of a decrease in productivity at the lowest altitudes and of an increase at the highest altitudes (as described by e.g., [51] for crop species; [21,22]). Accordingly, one may investigate future potential for pasture colonization at the highest altitudes, as limited by area availability and topography (e.g., [51]).
Our results are further coherent with former studies about suitable areas for pasture growth. With the increasing summer temperature, pasture species may find optimal temperature conditions at higher elevations [23]. Indeed, one notices that below 1100 m asl ( Figure 6), pastures may not find suitable conditions, and productivity may decrease. As reported in other cases (see e.g., [89,90]), Trisetum flavescens shows a potential decrease in presence and biomass in response to high temperatures in the summer, perhaps combined with a decrease in precipitation [90,91]. At higher elevations, above 2000 m asl, the large increase in productivity under climate change scenarios of the Nardus stricta may be explained by considering warmer temperature conditions and the properties of this species, which can adapt itself to drought conditions and needs less water than other pasture species [23,92].
Other studies demonstrated a potential increase in pasture productivity due to an increase in temperature and the subsequent reduction of snow cover duration. While in the CR period the growing season mostly overlaps with summer due to the optimal temperature range for growth, in the future potential anticipation of the GS, due to lack of snow cover, a sooner accumulation of thermal time, increasing densification and biomass of plants, anticipation of the peak of growth, and thereby an increase in productivity is seen [93,94].
From our calculated agroclimatic indices, (change in) annual pasture species productivity is driven by an increase in T during GS and less by P GS , which would not be a limiting factor in spite of the reduction of precipitation in some scenarios. As reported in other studies (e.g., [93,95,96]) grassland growth in Alpine catchments does not normally display water limitation, while low temperature at high altitudes may limit growth. In the results of this work, in future climate scenarios such limitation would be compensated by warming in all seasons. Water limitations may occur in some cases at the lowest altitudes, as driven by largely increased ET demand under temperature increase. Looking at Figure 7f,g, a small increase in ET efficiency and an increase in relative ET at low altitudes can be expected due to the reduction of precipitation. This explains the reduction of biomass for Trisetum flavescens, sensitive to drought periods.
An important novelty of this study is the use of updated IPCC AR6 scenarios and the comparison between AR5 and AR6.
In particular, we highlighted the larger increase in productivity and the more significant reduction of the specific water footprint for the AR6 scenarios.

Choice of Two Index Pasture Species
Clearly, in spite of the acceptable results here, improvements need to be pursued henceforth. The choice we made here of an index species (in each altitude zone) is clearly a simplified one, also taken in other studies [83,[85][86][87]. Two of the most abundant species in the area were chosen. As reported, species abundance was verified using some data available in fulfilment of the project MOUPA (Confalonieri R., Movedi E., personal communication, 2020), namely samples taken during summer 2019 in the area of Dosdé catchment (Northwestern Valtellina/Valdidentro) above 2000 m asl. Thereby, Nardus stricta was the first or second most abundant species. Moreover, in the experimental site of Alpe Boron, Nardus stricta was the second most abundant species. Looking at other studies about the abundance of pasture species, Nardus stricta is often found in Alpine high-altitude pastures [97,98]. This is related to the ability of the species to adapt to different climate and soil conditions [23]. At low altitudes, under 2000 m asl, Trisetum flavescens is often present in alpine pastures and managed permanent grassland, where cuts are applied during the growing season [53,54,99,100]. For this reason, this species was used here as an indicator of low-altitude pastures. In dealing with a single species, simulation is faster, and it is possible to consider a large area, such as the Sondrio district here.
We maintained the spatial support for simulation (i.e., the pasture areas in the valley) unchanged in the future. In doing so, we could analyze the productivity of our target species on the whole territory, in contrast to other studies focused on single farm system [87,101], and identify better or worse areas for growth. Moreover, it was possible to evaluate the behavior of single pasture species and their reaction to the changing climate to consider its potential for livestock.
Clearly, changes of pasture composition in response to competition between species, and alterations in response to climate change, including shifting to higher/lower altitudes and abandonment of some areas, may happen in the future.

Limitations and Future Improvements of the Simulation
Interspecific competition between different species competing for nutrients, water, and light and with different periods of growth may modify the dynamics of pasture, also in response to climate variability (e.g., [25]). The development of a version of the Poli-Pasture model able to simulate competition between multiple species is ongoing; however, it requires a large amount of information of pasture composition and dynamics and further refinement. Notice further that the dynamics of pasture cover in space may also depend upon anthropic management, and accordingly, changes of the spatial extent of simulation in the future may need some hypothesis about land cover evolution. Such topics clearly indicate lines worthy of investigation henceforward.
The Poli-Pasture model does not explicitly account for the direct impact of extreme rainfall (storm) events, if not indirectly via increased runoff, and work in this direction may be carried out. Further developments that could be discussed include assessment of the nutrient budget, given that the model assumes now full availability.

Conclusions
We analyzed the potential impacts of prospective climate change as projected under the IPCC AR5/AR6 scenarios upon two pasture and grassland species growth in the province of Sondrio in the central Italian Alps. We defined scores, or agroclimatic indices AIs, to objectively evaluate (i) climate suitability, (ii) pasture species productivity, and (iii) water use in the area.
An original contribution here concerns the use of a semi-distributed model for the simulation of high-altitude pasture areas over a large territory. Moreover, we projected future conditions considering new AR6 scenarios. Finally, in using AIs we could identify climate and hydrological conditions that could improve or worsen the growth of pasture species.
We conclude that an increase in mean temperature as expected during the 21st century may bring better conditions to the chosen pasture species on average by (i) increasing growing season length and (ii) increasing biomass. Uncertainty in future precipitation is seemingly less relevant for the growth in this water-abundant area, with likely no need of irrigation. An increased frequency of heat waves may occur, and a slight decrease in Trisetum flavescens biomass may be seen at the lowest altitudes (below 1100 m asl); however, this will be largely offset by an increased Nardus stricta yield at the highest altitudes.
Our results are of interest to scientists in the field of pasture/crop sciences willing to deepen into the potential impacts of climate change and to policy makers willing to explore strategies for adaptation. We further provided here possible hints for adaptation to breeders/farmers in Valtellina henceforth. Namely, one may suggest uplifting of pasturerelated activity at higher altitudes than now, to concentrate investments in an area where potentially larger productivity will be attained, possibly aiding economic development in the valley.