Exploring Short-Term Climate Change Effects on Rangelands and Broad-Leaved Forests by Free Satellite Data in Aosta Valley (Northwest Italy)

: Satellite remote sensing is a power tool for the long-term monitoring of vegetation. This work, with reference to a regional case study, investigates remote sensing potentialities for describing the annual phenology of rangelands and broad-leaved forests at the landscape level with the aim of detecting eventual effects of climate change in the Alpine region of the Aosta Valley (Northwest (NW) Italy). A ﬁrst analysis was aimed at estimating phenological metrics (PMs) from satellite images time series and testing the presence of trends along time. A further investigation concerned evapotranspiration from vegetation (ET) and its variation along the years. Additionally, in both the cases the following meteorological patterns were considered: air temperature anomalies, precipitation trends and the timing of yearly seasonal snow melt. The analysis was based on the time series (TS) of different MODIS collections datasets together with Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) collection obtained through Google Earth Engine. Ground weather stations data from the Centro Funzionale VdA ranging from 2000 to 2019 were used. In particular, the MOD13Q1 v.6, MOD16A2 and MOD10A1 v.6 collections were used to derive PMs, ET and snow cover maps. The SRTM (shuttle radar topography mission) DTM (digital terrain model) was also used to describe local topography while the Coordination of Information on the Environment (CORINE) land cover map was adopted to investigate land use classes. Averagely in the area, rangelands and broad-leaved forests showed that the length of season is getting longer, with a general advance of the SOS (start of the season) and a delay in the EOS (end of the season). With reference to ET, signiﬁcant increasing trends were generally observed. The water requirement from vegetation appeared to have averagely risen about 0.05 Kg · m − 2 (about 0.5%) per year in the period 2000–2019, for a total increase of about 1 Kg · m − 2 in 20 years (corresponding to a percentage difference in water requirement from vegetation of about 8%). This aspect can be particularly relevant in the bottom of the central valley, where the precipitations have shown a statistically signiﬁcant decreasing trend in the period 2000–2019 (conversely, no signiﬁcant variation was found in the whole territory). Additionally, the snowpack timing persistence showed a general reduction trend. PMs and ET and air temperature anomalies, as well as snow cover melting, proved to have signiﬁcantly changed their values in the last 20 years, with a continuous progressive trend. The results encourage the adoption of remote sensing to monitor climate change effects on alpine vegetation, with particular focus on the relationship between phenology and other abiotic factors permitting an effective technological transfer.


Introduction
Satellite-based remote sensing has been proved to effectively support vegetation monitoring [1]. Nowadays, satellite data are widely used, at the global scale, for bioclimatologic and bio-meteorologic issues, often related to the global climate change.
An ecological and evolutionary dilemma is posed to a variety of organisms, in particular vegetation, nowadays because of the environmental modifications related to climate change: can they respond properly to these ongoing changes? Does their phenology remain synchronous with other mechanisms or species that they interact with? Can they adapt their responses consistently with regard to the timing of reproduction? These questions are difficult to be answered without long-term observations and experiments. Natural/seminatural vegetation, such as rangelands and broad-leaved forests in mountain regions, can play an important role in exploring this topic. In fact, unmanaged vegetation, located in more extreme climatic situations (such as mountain rangelands and forests), is the first one that could be stressed, if a significant change affects local climatic conditions. Consequently, the long-term monitoring of rangelands behavior could drive one to better understand the dynamics of reaction/adaptation (resilience) of vegetation to climate change. In this context, remotely sensed data from EO (earth observation) missions operating over a longtime range could proficiently support this analysis. In the last few years, many works have been published in the literature presenting approaches, or case studies, about this topic with a significant improvement [2][3][4]. With special concern about mountain rangelands, their phenology at high altitudes is mainly related to snowpack that consequently, because of climate change, is drastically reducing [5]. Snowpack reduction depends on many factors related to the global, regional, and local climate. Global factors include changes of temperature and precipitation regimes-in alpine areas, the average temperature and rainfalls are increasing-moreover, rainfalls are tending to move from snowy to watery ones [6][7][8]. These environmental changes appear to have anticipated snowmelt, making spring begin earlier, as already observed in the Rocky Mountains (United States) [9] and in other mountain areas worldwide [10]. Moreover, phenological changes at high latitudes prove to be similar to the ones observed at high altitude [11]. The earlier beginning of the growing season, related to the anticipated snowmelt, can have multiple consequences: it could increase the length of the photosynthetic period, if the end of the season does not move or move to a later date. This impacts drought-related dynamics since an earlier snowmelt, and the consequent longer snow-free period, may increase the exposure of plants to this stress [12]. Additionally, an earlier snowmelt can significantly alter the dates on which species may come into bloom throughout the summer [13], because the ground and air will warm up earlier when snow disappears. For some species, a correlation was found between the timing of snowmelt and the abundance of flowering. In fact, one of the factors linking the dates of snowmelt to flowering abundance is frost [14]. If the probability of spring frost at a certain date remains the same and leaf or flower buds have anticipated, they become more vulnerable to frost: an increase in the frequency of frost-related damages is consequently expected for frost-sensitive species. Frost damage might also increase if the date of the last spring frost posticipates and if the rate of change in frost dates is slower than the one of snowmelt dates. Under these conditions, the knowledge of phenological trends at a wide scale within a certain area is crucial to adapt and plan pasture management under climate change. Land surface phenology (LSP) observations from remotely sensed data can provide key information related to seasonal dynamics [15] in the net ecosystem exchange of carbon [16] and surface energy and water balance [17] at regional-to-global scales [18]. LSP has also been shown to reflect the year-to-year variation of weather and is, therefore, a sensitive indicator of the ecosystem's response to long-term changes in the climate [19]. Consistently with perturbations expected from climate change [20], long-term phenology studies have identified widespread trends towards earlier springs and a lengthening of the growing season [21][22][23]. In contrast, observed trends of leaf senescence delay by remote sensing appearing less pronounced and being spatially heterogeneous are also supported [24,25]. All of these changes are expected to have important implications for ecosystem productivity and function [26].
Higher air temperatures have resulted in greater fluxes of precipitation (P) and evapotranspiration (ET) between terrestrial and atmospheric systems, leading to the intensification of the water cycle [27,28]. The growing season length partially controls ET, which cycles up to 62% of terrestrial water to the atmosphere [29]. P and the atmospheric evaporative demand have increased due to warmer temperatures [30], leading to changes in annual water yields [31], with humid regions tending to become wetter while arid regions have tended to become drier [32]. A longer growing season has the potential to alter water cycle fluxes through increased plant water use and ET [33] in the terrestrial ecosystems. The intensification of the water cycle is theoretically consistent with the Clausius-Clapeyron relation [34], which states that warmer air holds more water (for example, greater maximum specific humidity) and that consequently, atmospheric water vapor tends to increase non-linearly with increases in the air temperature. A warmer climate can also enlarge the growing season and provide more energy for ET. Therefore, accurate information regarding the start and the end of the vegetation growing season is important for monitoring and modeling the impacts of climate change on ecosystem processes. Studies exploring growing season variations can be effectively supported by both field observations and satellite data. A large number of ground-based studies have explored the phenological response of vegetation to climate change at particular locations or for subsets of species [35]. Satellitebased monitoring of phenology, on the other hand, allows the regional-to-global analysis of ecosystem responses to climate change, but shows some limitations while modeling and analyzing vegetation at the species level. Early remote sensing-based studies of phenology used coarse spatial resolution imagery from the advanced very high-resolution radiometer (AVHHR). In recent years, however, most of the studies were carried out at a higher spatial resolution (250, 500 m) with reference to image time series derived from the NASA's moderate resolution imaging spectroradiometer (MODIS) [36]. Most recently, a number of studies have begun to use 10-30 m resolution imagery, having similar spectral properties, available from Landsat and Sentinel 2 in order to improve phenological deductions at a higher spatial resolution [37][38][39][40][41][42][43][44].
A variety of different methods have been used to extract the so-called phenological transition dates, i.e., metrics defining the day of year associated with specific pheno-phase transitions [37], from time series of spectral vegetation indices such as the normalized difference vegetation index (NDVI) [39,40]. In particular, the start of season (SOS) and end of season (EOS) dates can be estimated from remotely sensed time series using a variety of methods to model observed temporal dynamics, including harmonic analysis [41], wavelet decomposition [42], splines [43], and curve fitting [45]. Curve fitting algorithms are the most commonly used, and among them, the asymmetric Gaussian, double logistic, and piecewise logistic functions are widely applied [46]. Regardless of the algorithm used, the overall methodology for extracting transition dates generally includes the following steps: satellite vegetation index time series generation; time series pre-process including gap filling, cloud screening, smoothing, the identification of valid phenological cycles, function fitting (usually splines or double logistic), and lastly, either parameter inversion or threshold date extraction [43,47].

Goals and Summary Description
Under these premises, this work was aimed at exploring short-term climate change effects, by remotely sensed data, on phenological metrics (PMs) and evapotranspiration (ET) of rangelands and broad-leaved forests in Aosta Valley, an alpine region located in north-western Italy. In order to analyze possible climate change effects in this region, particular attention was paid to natural and semi-natural vegetation classes: rangelands (R) and broad-leaved forests (B). The expectation was that climatic factors affecting vegetation phenology were snowpack time persistence (snow cover duration), precipitation and air temperature anomaly. Specifically referring to alpine rangelands and broad-leaved forests, only a few papers in the literature [10] try to relate the effects of climate change on vegetation PMs and ET with snowpack time duration, air temperature anomalies and precipitation. This is probably due to the fact that ground measures in mountain regions are difficult to be obtained with a proper spatial distribution. Moreover, only in recent years scientific community attention has been focused on mountains, thanks to their role of climate change sentinels. This work, based on freely accessible remotely sensed data, is intended to fill this gap with special concerns about the Alps.
It is worth remembering that rangelands have a considerable importance for both ecosystem and economic services of this region, since livestock is one of the main sources of economic income. Moreover, rangelands have a wide spatial distribution over the territory; therefore, monitoring their phenology and physiology can certainly support the comprehension of the local climate and territorial dynamics. Differently, broad-leaved forests are known to be particularly sensitive to climate change, being located in two very sensitive zones: the bottom of the valley (where precipitations and temperature are highly varying) and the higher altitude bands (for e.g., Alnus alnobetula (Ehrh.) K. Koch and Acer pseudoplatanus L.; Fraxinus excelsior L. etc.), where variations in snow cover time, persistence and temperature can deeply affect vegetation phenology.
Additionally, together with conifers (and in particular Larix decidua Mill.) they are greatly interesting for the National Inventory of Forests and Forest Carbon pools (INFC), and, at the regional level, for forestry planning adaptation to the climate change perspective.
With these premises and with reference to the CORINE land cover map (see forward on), the analysis was focused on pastures, natural grasslands, moors/heathlands and broad-leaved forests. To explore the altitude effects on trends, the analysis was achieved at different height class. PMs/ET investigations were based on TERRA MODIS image time series (MOD13Q1 and MOD16A2 products, see forward on) ranging from 2000 to 2019. With reference to the MOD13Q1 product, a NDVI image time series was generated and used to estimate the PMs: SOS (start of season), EOS (end of seasons), LOS (length of season), MAXVI (maximum vegetation index). In order to validate the results concerning PMs' estimates as computed by the TIMESAT tool, the MCD12Q2 phenological collection was adopted.
Air temperature anomalies, snowpack timing and precipitation data were obtained from both EO datasets and ground stations.
With reference to the MOD16A2 product, the corresponding ET (computed as yearly cumulated eight days observation per year) image time series was generated and used to describe trends affecting rangelands and broad-leaved forests at different altitudes. The results showed that most of the above-mentioned metrics significantly changed in the last 20 years in this region, proving that this alpine area is very sensitive to climate change effects. The trends were finally modelled by first order polynomial relationships and some operational information was extracted and discussed, with the aim of providing a forecasting model useful for rangeland management under the climate change conditions.

Study Area
The study area corresponds to the entire Aosta Valley, an administrative region located in the north-western Alps in Italy, close to France and Switzerland. It is an alpine region ( Figure 1 Table 1 reports the altitude ranges and their spatial distribution, respectively. South and the Grand Combin (in Switzerland) in the North. Table 1 reports the altitude ranges and their spatial distribution, respectively.
Given its topographic features, mostly related to the local presence of the highest European mountains, it can be reasonably considered an open-cell laboratory or environmental and biological sciences.  The Aosta Valley climate is strongly influenced by topography, as well. Surrounding high mountains prevents the access of humid air masses from Mediterranean or Atlantic areas, making local climate characterized by: a high degree of aridity in the central area (rainfalls lower than 500 mm·y −1 ); an average rainfall (above 1400 mm·y −1 ) in the border  Given its topographic features, mostly related to the local presence of the highest European mountains, it can be reasonably considered an open-cell laboratory or environmental and biological sciences.
The Aosta Valley climate is strongly influenced by topography, as well. Surrounding high mountains prevents the access of humid air masses from Mediterranean or Atlantic areas, making local climate characterized by: a high degree of aridity in the central area (rainfalls lower than 500 mm·y −1 ); an average rainfall (above 1400 mm·y −1 ) in the border areas and especially in its south-eastern and north-western valleys [48]. In winter, precipitations are mainly snowy. In summer, convective rainfalls are quite common, determining frequent thunderstorms; spring and autumn are characterized by stratified rainfalls, possibly lasting for several days, with a consequent increase in the flood risk. In force of its highly variegated climatic situation, Aosta Valley land cover is heterogenous.

Remote Sensing Data
MOD13Q1 v.6 [49] product from the NASA TERRA MODIS (moderate resolution imaging spectroradiometer) mission was used to map the normalized difference vegetation index (NDVI) over the area in the period from 18 February 2000 to 31 December 2019. Four hundred-fifty-two MOD13Q1 images were obtained for free by Google Earth Engine [50]. They were stacked into a NDVI time series (NTS) and filtered according to the pixel reliability layer, supplied together with the product, to map reliability of each scene pixel. Stacking and filtering (Savitzky-Golay, [51,52]) were operated by a self-developed IDL 8.1 routine [53]. The vegetation index from the MOD13Q1 v.6 product is a composite one, having a time step of 16 days and a spatial resolution of 250 m. The MOD13Q1 v.6 product provided a vegetation index (VI) value at a per pixel basis. It is referred to as a continuity index to the existing one from the National Oceanic and Atmospheric Administration advanced very high resolution radiometer (NOAA-AVHRR) data. A composition algorithm selected the best available pixel value from all the acquisitions within the considered 16-day period. The selection criteria are the following: no clouds, low viewing angle, and highest NDVI value, adopting the classic formula as follows: (Equation (1)).
where ρNIR and ρRED are the at-the-ground reflectance of the near-infrared (MODIS band 2) and red (MODIS band 1) bands, respectively. ET (evapotranspiration) maps were retrieved from the MOD16A2 v.6 collection [54] by Google Earth Engine for the period from 1 January 2000 to 31 December 2019 for a total of 868 images.
The MOD16A2 v.6 Evapotranspiration/Latent Heat Flux product is an 8-day composite product with a geometric resolution of 500 m. The algorithm adopted for the MOD16A2 data product collection refers to the Penman-Monteith equation [55], which includes inputs of daily meteorological reanalysis data such as albedo, land cover, vegetation property dynamics and LST (land surface temperature) that in some studies were adopted to evaluate ET [56] or heat fluxes [57]. Pixel values for the net evapotranspiration (ET) is the sum of all 8 days within the composite period expressed in kg m −2 8d −1 . The MOD16A2 v.6 ET layers were stacked into an ET time series (ETS); no filtering was applied. Stacking was operated by a self-developed IDL 8.1 routine.
The Shuttle Radar Topography Mission (SRTM) digital elevation model was used for this work [58,59]. The SRTM v3 product (SRTM Plus) was provided by NASA JPL with a grid size of 1 arc-second (approximately 30 m). Native SRTM was resampled at 30 m (SRTM30) to make pixel squared.
The MOD13Q1 v.6 and MOD16A2 v.6 layers were geometrically oversampled up to 30 m to refine the area zonation by applying the bicubic spline interpolation [60] method in order to have grids with the same spatial resolution. In the case of a lack of data, an interpolation was performed in R v.1.3.1093 by adopting the package "Impute TS" specializes on univariate time series imputation. It offers multiple state-of-the-art imputation algorithm implementations along with plotting functions for time series missing data statistics. The following functions were adopted in R Studio: ggplot_na_distribution is a function that visualizes the distribution of missing values within a time series. The time series is plotted and whenever a value is NA, the background is colored differently. This gives a nice overview as to where in the time series, most of the missing values occur.
ggplot_na_imputations is a function that plots the imputed values for a time series. Therefore, the imputed values (filled NA gaps) are shown in a different color than the other values [61].
To validate PMs estimates from the NDVI of MOD13Q1 v.6 collection, the MCD12Q2 v.6 Land Cover Dynamics product (informally named the MODIS Global Vegetation Phenology product) was adopted. It was obtained from Google Earth Engine with a 500 m spatial resolution covering the period from 1 January 2001 to 31 December 2018. Missing values were interpolated, or extrapolated, for a maximum of one year before or after the reference time period of the collection by adopting the functions ee.interpolate and ee.Model.fromAiPlatformPredictor available in Google Earth Engine libraries. The dataset provides estimates of the timing of vegetation phenology at global scales. Additionally, it provides information related to the range and summation of the enhanced vegetation index (EVI) computed from MODIS surface reflectance data at each pixel. It identifies the onset of greenness, green-up midpoint, maturity, peak greenness, senescence, green-down midpoint, dormancy, 2-band enhanced vegetation index (EVI2) minimum, EVI2 amplitude, integrated EVI2 over a vegetation cycle, as well as overall and phenology metric-specific quality information. The MCD12Q2 v.6 data product is derived from time series of the 2-band enhanced vegetation index (EVI2) calculated from the MODIS nadir bidirectional reflectance distribution function (BRDF)-adjusted reflectance (NBAR). Vegetation phenology metrics are identified for up to two detected growing cycles per year. For pixels with more than two valid vegetation cycles, the data represent the two cycles with the largest NBAR-EVI2 amplitudes [62].
All data were converted from native geographical reference systems WGS84 into the ED50 UTM 32N one.

CHIRPS (Climate Hazards Group InfraRed Precipitation with Station data) is a 30+
year quasi-global rainfall dataset. CHIRPS incorporates 0.05 • (5550 m) resolution satellite imagery with in-situ station data to create gridded rainfall time series for trend analysis and seasonal drought monitoring [63]. A time series of precipitation estimates was obtained for the study area by Google Earth Engine. The annual rainfall of the whole Aosta Valley Region was considered along time. No spatial resampling was performed, and the native spatial resolution maintained.
NDSI Snow Cover: the MOD10A1 v.6 Snow Cover Daily Global 500 m product contains snow cover, snow albedo, fractional snow cover, and quality assessment (QA) layers. Snow cover data are based on snow mapping algorithm that employs a normalized difference snow index (NDSI) and other criteria tests [64].

Land Cover Data
To properly describe land cover at regional level, the CORINE land cover 2018 dataset (hereinafter referred to as CLC2018) was used; the CORINE (Coordination of Information on the Environment) land cover inventory was initiated in 1985 to standardize data collection on land in Europe to support environmental policy development. The project is coordinated by the European Environment Agency (EEA) in the frame of the EU Copernicus program and implemented by national teams. The number of participating countries has increased over time, currently including 33 member countries and six cooperating countries with a total area of over 5.8 million km 2 . CLC2018, specifically, is one of the available datasets produced within the general CORINE frame and refers to land cover/land use status in 2018. The reference year of the first CLC inventory was 1990 and the first update was achieved in 2000. The current update cycle is 6 years. Satellite imagery provides the geometrical and thematic basis of maps; it is integrated with in-situ data as essential ancillary information. The CLC Level 3 dataset considers 44 classes with a minimum mapping unit (MMU) of 25 hectares (100 m minimum mapping width). The CLC2018 Level 3 of Aosta Valley was obtained from the collection COPERNICUS/CORINE/V20/100 m by Google Earth Engine in grid format (100 m grid size) [65][66][67]. The CLC2018 Level 3 classes and corresponding statistics are reported, for the whole Aosta Valley, in Table 2.
Land cover classes were considered with reference to the considered 3 height classes (see forward on), in order to give a more comprehensive description on the ongoing processes affecting vegetation and the role of altitude in vegetation responses. As indicated in Table 2 (green cells), this study focused on broad-leaved forests and rangelands, that ware assumed to correspond to the following CORINE land cover classes: pastures, broad-leaved forests, natural grasslands, moors and heathlands ( Figure 2). B and R were chosen because they both are very sensitive to climate change and have huge importance in Aosta Valley in terms of the economic benefits and the great biodiversity hosted.

Ground Data (Weather Stations Data)
In order to analyze the possible relationship of PMs and ET with climate change, the air temperature and precipitation data were obtained for the explored period (2000-2019) from the Data Center of the Centro Funzionale della Valle d'Aosta (Functional Center of Aosta Valley). The dataset includes temperature values as collected by 100 weather stations distributed all around the Aosta Valley region [68]. For each station, the daily mean air temperature was considered. Stations were preliminarily assigned to the three altitude classes (namely H1, H2, H3 see forward on) that authors considered for this work. Only stations falling in rangelands and in proximity of broad-leaved forests or sufficiently close to them (distance < 3 km) were considered. Furtherly, they were reduced to guarantee that  Aosta Valley). The dataset includes temperature values as collected by 100 weather stations distributed all around the Aosta Valley region [68]. For each station, the daily mean air temperature was considered. Stations were preliminarily assigned to the three altitude classes (namely H1, H2, H3 see forward on) that authors considered for this work. Only stations falling in rangelands and in proximity of broad-leaved forests or sufficiently close to them (distance < 3 km) were considered. Furtherly, they were reduced to guarantee that the same number of stations (finally, 15) were assigned to the three altitude classes for a total of 43.

Methodology
The analysis was aimed at testing if any significant climatic trend affected phenological metrics (PMs, see forward on) and ET (as measured by NTS and ETS, respectively) in the period from 2000-2019 and if the time series (TS) observed were related to rising temperature, precipitation trend and the timing of annual seasonal snowmelt. CLC2018 classes (231 = pastures, 311 = broad-leaved forests, 321 = natural grasslands, 322 = moors/heathlands) were considered and aggregated into a single macro-class, hereinafter named "H".

Testing PM and ET Trends from NTS and ETS
To investigate changes affecting the vegetation activity in consequence of cli change, some metrics were extracted from the above-mentioned time series (NTS ETS). Phenological metrics (PMs) are numerical indices of vegetation activity, corresp ing to some key points identified along its annual growing season. Climate change pr to condition such activity shifting and reshaping of past "ordinary" behavior of p

Testing PM and ET Trends from NTS and ETS
To investigate changes affecting the vegetation activity in consequence of climate change, some metrics were extracted from the above-mentioned time series (NTS and ETS). Phenological metrics (PMs) are numerical indices of vegetation activity, corresponding to some key points identified along its annual growing season. Climate change proved to condition such activity shifting and reshaping of past "ordinary" behavior of plants along the year. The following PMs were considered: the start of the growing season (SOS) representing the day of the year (DOY) when phenology is admitted to boost; the end of the season representing the day of the year (DOY) when phenology is admitted to stop; the length of the growing season (LOS) representing the time range (in number of days) separating EOS from SOS; the maximum of NDVI (MAXVI) representing the highest value reached by NDVI during the growing season and proved to be a good predictor of climate change effects on vegetation.
PMs were estimated by TIMESAT 3.3 with STL software [72][73][74][75] that was specifically developed to enable the monitoring of land surface processes by remotely sensed data. TIMESAT 3.3 with STL [76][77][78] iteratively fits and smooths by mathematical functions the yearly NDVI time series, finding the best smoothed approximation of the NDVI along the year at the pixel level. Once raw data have been approximated by the selected fitting function, PMs can be extracted in correspondence to singular points having a phenological meaning (e.g., EOS, SOS, LOS, MAXVI, etc.) along the local temporal profile. In this work, phenological metrics (SOS, EOS, LOS, MAXVI) were estimated at the pixel level by the TIMESAT 3.3 software tool processing the whole NTS. A seasonal trend decomposition by Loess (STL) was adopted to de-trend the NTS pixel profile and remove noise. The seasonal component was refined by Savitzky-Golay filtering to reduce, but not remove, the remaining local strong variations. The yearly growing season was recognized with the whole multi annual time series using the sinusoidal harmonics approach.
With reference to de-trended/filtered NTS profiles, PMs were extracted using the simple relative thresholding approach. The approach assumes that the SOS and EOS correspond to the date (DOY) when the NDVI value reaches a quote (0.5 for this work) of the range [0.2-MAXVI] during the starting and ending phenological phases, respectively. The LOS map was computed by grid differencing from EOS and SOS maps. The MAXVI was found looking for the highest NDVI value (yearly averaged) observed between the SOS and EOS. In order to test the quality of the PMs obtained, they were validated by comparing those obtained from the MODIS phenology product MCD12Q2, with an overall discrepancy between them being less than 5%.
With reference to ETS and separately for H1, H2, H3, the class yearly average ET value was computed and plotted along the years (2000-2019) at regional level, to look for trends.
The analysis was performed at the pixel and year level. Consequently, estimates of all PMs were represented as raster layers, mapping their spatial distribution at the considered year.
Considering both PMs and ET estimates, a first investigation was aimed at removing outliers from data. This was obtained by calibrating a first order polynomial model relating PMs/ET with the years and testing the residuals (model predicted-observed).
where M is the generic PM/ET yearly averaged metric and Y the considered year (in this case of study, the period between 2000-2019).
Data showing an absolute percent residual (ε, Equation (3)) > 100% were labelled as outliers and removed. A new model calibration was consequently run after removal.
Scatterplots are reported in Figure 4; model coefficients and coefficient of determination (R 2 ) and p-values (calculated by PAST software [69]) are reported in Table 4A,B. In order to test the quality of the PMs calculated from the NDVI of MOD13Q1 v.6 collection, the MCD12Q2 v.6 Land Cover Dynamics product (informally known as the MODIS Global Vegetation Phenology product) was adopted to perform the validation.   Reported values correspond to the yearly mean value of the 8 days aggregated data from MOD16A2 product. In all the graphs orange points represent outliers (see Equation (2)). All the trends proved to be significant (see Table 4B). DOY: day of the year.
According to Equation (3), the number of outliers was found to be lower than four (orange points in Figure 4) for all considered parameters (PMs and ET). After, the outlier removal of PMs and ET trends turned to be statistically significant.  (H1, H2, H3) yearly mean values of ET. Reported values correspond to the yearly mean value of the 8 days aggregated data from MOD16A2 product. In all the graphs orange points represent outliers (see Equation (2)). All the trends proved to be significant (see Table 4B). DOY: day of the year.

Yearly Snowmelt-Snow Cover Timing
The timing of the annual seasonal snowmelt and any related potential change has broad ecological implications and can impact vegetation and human livelihoods, particularly in and around high latitude areas and mountainous systems. For this task, the code available in Google Earth Engine (https://code.earthengine.google.com/2db2624ad1 0daba32e1155e0dfc4e03a (accessed on 26 January 2021)) was used. It made it possible to produce annual maps reporting the first day within a year where a given pixel reaches zero percent of snow cover and the corresponding snow cover timing.

Meteorological Data Processing: Ground Stations and CHIRPS
Available daily long-term meteorological data from ground stations of the air temperature were preventively homogenized, to ensure that their variability was due to climatic factors, solely. HOMER and Splidhom methods were used [79,80] and the homogenized series were tested through a quality control process aimed at detecting anomalous records [81]. This process was performed in RStudio v.1.3.1093. Finally, the mean daily air temperature was averaged along the year separately for the three considered altitude classes. Additionally, the CHIRPS data were geometrically oversampled up to 30 m to refine area zonation by applying the bicubic spline interpolation [60] method in order to have grids with the same spatial resolution of the entire collection adopted. The CHIRPS were used only to observe how precipitations were temporally and spatially distributed in the whole territory over the time. For the specific analysis presented in this work, considering altitude bands and land cover (B and R), ground weather stations were adopted because they were spatially more suitable to the different dataset used and the kind of analysis performed.

Results
PMs were mapped over the area as raster layers and spatially averaged with respect to H1, H2 and H3 classes (Table 3) by ordinary zonal statistics.
Graphs of Figure 4a-c show that the average PMs values of H1, H2, and H3 significantly changed their values in the last 20 years, with a continuous progressive trend observable for all of them. First order polynomials used to model trends were calibrated after outlier removal. The coefficients are reported in Table 4, together with the corresponding coefficients of determination (R 2 ) and p-value. To operationally synthesize information about strength of changes, in Table 5, yearly rates of increment/decrement of PMs/ET are reported together with total increment/decrement rate in the period from 2000-2019, as derived from the calibrated models. Table 5. Yearly average and cumulated (2000-2019) variations of PMs and ET along the years as estimated by 1st order polynomial regression (Figure 4). Positive values indicate a delay or lengthening of season, negative values an advance of season according to the observed trends reported in Figure 4. The first part of the table indicates the yearly gain, while the second part (cumulated) represent the overall changes in the observed period analyzed for each PMs from the starting in the year 2000 up to the end (2019). Standard deviations have been calculated to define the uncertainty margin on the basis of the estimated gain (indicated after the ±). According to Equation (3), the number of outliers was found to be lower than four (orange points in Figure 4) for all considered parameters (PMs and ET). After, the outlier removal of PMs and ET trends turned to be statistically significant.

RANGELANDS
The following interpretation of results can be given. In rangelands, the start of the growing season (SOS) significantly anticipated every year during the entire considered period (2000-2019). The yearly anticipation of SOS was estimated to be about 2-3 days per year depending on the altitude class (Table 5). Conversely, the end of the season (EOS), significantly delayed every year during the entire considered period (2000-2019). The yearly delay of EOS was estimated to be about 2.5-3 days per year depending on the altitude class (Table 5). Consequently, the length of the growing season (LOS) significantly increased every year during the entire considered period (2000-2019). The yearly enlargement of LOS was estimated to be about 5-6 days per year depending on the altitude class (Table 5). Similarly, MAXVI showed an increasing trend too, but the yearly difference value was found to be not significant (see Discussion section). In broad-leaved forests, it is possible to observe that the start of the growing season (SOS) significantly anticipated every year during the entire considered period (2000-2019). The annual beginning of the season (SOS) was estimated to be about 1 day per year considering all the altitude classes (Table 5). Conversely, the end of the vegetative season (EOS), significantly delayed every year during the entire considered period (2000-2019). The yearly delay of EOS was estimated to be about 2 days per year in all altitude classes (Table 5). Therefore, the length of the growing season (LOS) significantly increased every year during the entire considered period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). The yearly enlargement of LOS was estimated to be about 3.5 days per year per each altitude class (Table 5). Similarly, MAXVI showed an increasing trend too.
The obtained results showed that all PMs and ET significantly changed in the reference period for both broad-leaved forests and rangelands in Aosta Valley. In particular, a progressive linear enlargement of rangelands and softwood growing season was observed for all the considered altitude classes (H1, H2, H3).
Evapotranspiration (ET) significantly increased every year during the entire considered period (2000-2019) for both rangelands and broad-lived forests. In rangelands, the yearly increase in ET was estimated as 0.06 Kg m −2 (H1), 0.04 Kg m −2 (H2) and 0.06 Kg m −2 (H3) per year, respectively ( Table 5). While in broad-leaved forests, the annual increase in ET was estimated as 0.08 Kg m −2 (H1), 0.05 Kg m −2 (H2) and 0.06 Kg m −2 (H3) per year, respectively ( Table 5). In Table 5, the cumulated differences computed for the whole reference period are also reported.
In order to better understand what factors had strongly influenced the phenology, not only temperature and evapotranspiration were considered. In particular, in this analysis, precipitation and the snowpack persistence were also considered by defining the timing of the annual seasonal snowmelt; that is to say, the yearly first day of no snow cover ( Figure 5). By observing the precipitation (including rainfalls and snowfalls) trends in the whole Aosta Valley autonomous region, in the reference period of study, no significant variation or particular trend was observed with a stationary mean trend; nevertheless, a statistically significant decreasing trend was only observed in the H1 class. It is worth remembering that the principal central valley was affected in general by less precipitation than the lateral valleys. In fact, the water supply comes climatically and historically from the side valleys, where rainfall is greatest. These evidences find testimony in the numerous hydraulic works of derivation towards the main valley [82]. This evidence is particularly important in regard to the present and future adaptation of the perspective of vegetation management to climate change. The obtained results showed that all PMs and ET significantly changed in the reference period for both broad-leaved forests and rangelands in Aosta Valley. In particular, a progressive linear enlargement of rangelands and softwood growing season was observed for all the considered altitude classes (H1, H2, H3).
Evapotranspiration (ET) significantly increased every year during the entire considered period (2000-2019) for both rangelands and broad-lived forests. In rangelands, the yearly increase in ET was estimated as 0.06 Kg m −2 (H1), 0.04 Kg m −2 (H2) and 0.06 Kg m −2 (H3) per year, respectively ( Table 5). While in broad-leaved forests, the annual increase in ET was estimated as 0.08 Kg m −2 (H1), 0.05 Kg m −2 (H2) and 0.06 Kg m −2 (H3) per year, respectively (Table 5). In Table 5, the cumulated differences computed for the whole reference period are also reported.
In order to better understand what factors had strongly influenced the phenology, not only temperature and evapotranspiration were considered. In particular, in this analysis, precipitation and the snowpack persistence were also considered by defining the timing of the annual seasonal snowmelt; that is to say, the yearly first day of no snow cover ( Figure 5). By observing the precipitation (including rainfalls and snowfalls) trends in the whole Aosta Valley autonomous region, in the reference period of study, no significant variation or particular trend was observed with a stationary mean trend; nevertheless, a statistically significant decreasing trend was only observed in the H1 class. It is worth remembering that the principal central valley was affected in general by less precipitation than the lateral valleys. In fact, the water supply comes climatically and historically from the side valleys, where rainfall is greatest. These evidences find testimony in the numerous hydraulic works of derivation towards the main valley [82]. This evidence is particularly important in regard to the present and future adaptation of the perspective of vegetation management to climate change. Finally, by analyzing snowpack persistence defined by the timing of the annual seasonal snowmelt that is calculated as the yearly first day of no snow cover (SC) in all the Finally, by analyzing snowpack persistence defined by the timing of the annual seasonal snowmelt that is calculated as the yearly first day of no snow cover (SC) in all the altitudes bands of the region, we observed a general anticipation trend in snowmelt. In all classes H1, H2, and H3, a linear regression was performed and all the trends were statistically significant, as shown by the p-value retrieved with the Mann-Kendall tests [83,84]. In particular, in H1 (in the bottom valley) the annual first day of no snow cover trend was particularly marked. A possible explanation can be found both in rising temperature anomalies and in the reduction in precipitation, as shown in the present study.
With respect to Aosta Valley, it was observed that between 2000 and 2019 the SOS was anticipated by almost 1-2 months in rangelands (R) and 20-25 days in broadleaved forests (B), the EOS was delayed by about 1.5 months for both B and R and the LOS stretched to about 3-4 months for all the considered altitude classes for both B and R. A cumulated significant increase was observed for both MAXVI and ET, as well. One of the possible factors that could be responsible for this important variation can be reasonably found in the increase in the air temperature that affected both snowpack duration (as described above) and vegetation phenology (whereas water is not a limiting factor). This was confirmed by the analysis of the annual anomaly of the average yearly air temperature in the period from 2000-2019. The air temperature values were obtained from the automatic weather stations and aggregated with respect to the considered altitude classes (H1, H2 and H3). The air temperature class anomalies showed an increasing and statistically significant trend, which increased with the altitude class ( Figure 6). Finally, a validation of the PMs retrieved on TIMESAT was performed. To do this, the MODIS product MCD12Q2 in Google Earth Engine was used. It is worth remembering that the MODIS product PMs are obtained starting from the EVI (enhanced vegetation index) and from the product GreenUp, Senescence, EVI_Amplitude and others, to consider SOS, EOS and LOS, respectively. Here, the relative trends between SOS, EOS and LOS are reported (Figure 7).

Discussions
PMs and ET modeling and trend analysis based on EO data, with particular concern about multispectral passive remote sensing, proved to provide useful information at seasonal, inter-annual and annual time scales. This can be proficiently used to describe phenomena related to biosphere-atmosphere exchanges of carbon, energy, and water at both regional and global scales. Mountain environments are known to be very sensitive to climate change, in particular to temperature increases [24]. Similar to other alpine areas, Aosta Valley is suffering from an important increase in temperature as a result of climate change, that appears to determine significant changes in the growing season of rangelands. The air temperature increase depends on altitude. The results showed that, in the analyzed period, the average increase in the yearly average temperature was lower than 2 °C, lower than 2.5 °C and higher than 3 °C within altitude classes H1, H2 and H3, respectively. These evidences can be directly related to the described trends of the considered phenological metrics and evapotranspiration. High altitudes are, in fact, more sensitive to rising temperatures that determine a reduced persistence of snowpack ( Figure 5) and a lengthening of the growing season, as pointed out in the analysis section. For this reason, they appear among the most severely and rapidly impacted ecosystems, being possibly impacted by any change in temperature and precipitation patterns at all scales [27]. Snow and ice are the main control parameters of the hydrological cycle, particularly of the seasonal runoff, and their variations can condition the entire geosystem (rocks, soils, vegetation, and river discharges). With climate change, water will probably become less available due to global warming and consequences will reach far beyond mountain regions [16]. Similarly, climate change is likely to increase exposure to either natural or economic hazards, all the more so because in many mountain areas, low incomes levels are higher than in lowland areas [7].

Discussion
PMs and ET modeling and trend analysis based on EO data, with particular concern about multispectral passive remote sensing, proved to provide useful information at seasonal, inter-annual and annual time scales. This can be proficiently used to describe phenomena related to biosphere-atmosphere exchanges of carbon, energy, and water at both regional and global scales. Mountain environments are known to be very sensitive to climate change, in particular to temperature increases [24]. Similar to other alpine areas, Aosta Valley is suffering from an important increase in temperature as a result of climate change, that appears to determine significant changes in the growing season of rangelands. The air temperature increase depends on altitude. The results showed that, in the analyzed period, the average increase in the yearly average temperature was lower than 2 • C, lower than 2.5 • C and higher than 3 • C within altitude classes H1, H2 and H3, respectively. These evidences can be directly related to the described trends of the considered phenological metrics and evapotranspiration. High altitudes are, in fact, more sensitive to rising temperatures that determine a reduced persistence of snowpack ( Figure 5) and a lengthening of the growing season, as pointed out in the analysis section. For this reason, they appear among the most severely and rapidly impacted ecosystems, being possibly impacted by any change in temperature and precipitation patterns at all scales [27]. Snow and ice are the main control parameters of the hydrological cycle, particularly of the seasonal runoff, and their variations can condition the entire geosystem (rocks, soils, vegetation, and river discharges). With climate change, water will probably become less available due to global warming and consequences will reach far beyond mountain regions [16]. Similarly, climate change is likely to increase exposure to either natural or economic hazards, all the more so because in many mountain areas, low incomes levels are higher than in lowland areas [7].
The results from this work demonstrate that phenological and evapotranspiration processes are significantly changing in the study area and that these effects are linked to rising temperatures and the timing of annual snowmelt. Changes can be modelled by linear trends, confirming what other similar works already found for different areas [12][13][14][15][16][17][18].
In particular, in Aosta Valley, PMs and ET and snow time were proved to have significantly changed their values in the last 20 years, with a continuous progressive trend observable for all of them.
As previously discussed in the analysis section, it is worth remembering that in terms of strength of changes, in rangelands, an average delay of the EOS was observed of about 2.6 days, while in broad-leaved forests of about 2.3 days, independent from the altitude class. Differently, the SOS in R appears to anticipate by about 2 d/y up to 2000 m ASL. (H1 and H2) and by about 3 d/y at higher altitudes (H3) while in B it appears to anticipate by about 1 d/y at all altitudes considered H1, H2, and H3. Consequently, in R, the LOS proved to enlarge by about 4.7 d/y at lower altitudes (H1 and H2) and by about 6.5 d/y at higher altitudes (H3) while in B proved to enlarge by about 3.8 d/y in all classes. According to Borgogno et al. [85], the potential accuracy of NDVI measurements is about ±0.02; consequently, the estimated yearly variations of MAXVI cannot be considered singularly significant. Nevertheless, the cumulated effects along the entire explored period (2000-2019) showed that MAXVI significantly changed, since the accuracy reference value of ±0.002 was largely overcome. MAXVI variations between 2000 and 2019 appear to be positively higher at lower altitudes (about +0.09) while almost stable as the altitude increases. This could be possibly explained with reference to the previously demonstrated increase in both biomass production (MAXVI) and enlargement of the growing season, that consequently make vegetation need more water yearly. This can be explained by admitting that at lower altitudes in Aosta Valley, grasslands and pastures are often irrigated [86,87]. Consequently, farmers can vary water release regimes to face climate change effects (higher temperatures, in particular) with the result of moving forage yields (that NDVI is a predictor of) to higher values. Differently, where more natural (not managed) systems are located (higher altitudes) the increase in yearly MAXVI can only be related to glacier melting that could compensate the increase in water requirement (as confirmed by the ET analysis) by vegetation: glaciers are, in fact, dramatically reducing in Aosta Valley and as shown here, snow cover time (annual snow melting) is also particularly evident and remarkable. Moreover, another compensating action could come from the surrounding forest areas that have been proved to tolerate summer heatwaves. These results find strong evidence in different studies in the literature [88][89][90][91][92][93][94][95][96][97][98][99]. Concerning on the variation between PMs in H1, H2 and H3 it can be said that in rangelands at different altitudes, the SOS had a more marked variation at high altitudes as can be seen from the gains (in Table 4A) as in the past the snow normally lasted longer, normally at the top rather than at the bottom of the valley. While in regard to EOS there is a more marked gain, as it is reasonable to expect at low altitudes than at high altitudes because of the lower temperatures, the frost and the strong temperature variations at altitude (and possible snowfalls) make the vegetative conditions less favorable. It should also be remembered that at low altitudes, the rangelands represent vegetation with a marked anthropic conditioning (such as the possibility of irrigation) which does not happen at high altitudes (except in H2 but in the height of summer) which makes vegetation a more exploitable resource, long given by the ease of access given to its topography and micro-climate (at the bottom of the valley). In regard to LOS, a more marked gain is observed at high altitudes due to the fact that the shorter duration of the snow and the increase in temperatures favor a more marked lengthening of the season at altitude than in the past (certainly due to the variation of the SOS and EOS). Finally, MAXVI has a more marked gain at low altitudes (H1), for the simple fact that the possibility of conditioning the vegetation with irrigation and agronomic practices significantly affects its productivity and relative vigor. In the case of broad-leaved forests in all classes, a similar gain can be observed in the SOS, EOS and LOS. This is due to the fact that trees are less sensitive and more adapted to environmental variation. Nevertheless, it is interesting to see that in H2 at middle altitude, the MAXVI gain is higher than in H1 and H3, probably because they found a more favorable ecological condition, but more studies can be carried out to find a good answer.
With reference to the ET graphs of Figure 4d,h, significant increasing trends were observed for all height classes. According to the values of Table 4, it can be noted that all height classes behave similarly; water requirement, in a period of 8 days, appears to averagely increase by about 0.05 Kg·m −2 (about 0.5%) every year, for a total increase of about 1 Kg·m −2 in 20 years (2000-2019), corresponding to a percentage difference in water requirement from vegetation of about 8%. This could possibly be explained by the increase in biomass production (well represented by MAXVI) and by the enlargement of the growing season.
The authors are conscious that this work just focuses on a semi-natural component, as rangelands and broad-leaved are. More investigations should be carried out with reference to other land cover classes and with more annual and interannual data at higher resolution. Nevertheless, preliminary results from this work could support breeders, foresters and policymakers to better address their alpine activities such as grazing and management of natural (water included) resources. A radical change is also expected by technicians and institutional subjects in their ordinary procedures, by adapting the scheduling of alpine activities according to this information. It is the authors' opinion that the greatest limitation to the technology transfer of this type of approach is a lack of "georeferenced" data from ground activities needed to calibrate and validate EO data-based deductions. A great effort is therefore required from involved players to carefully consider the possibility of structuring proper networks for data acquisition and transfer and procedures where georeferencing of ground observation is mandatory. Georeferencing of ground data is, in fact, at the basis of an effective and reliable spatial-based approach such as the one here proposed with a time series, where EO data (especially if available over a long-time span) play a crucial role. Consequently, we expect that future approaches should more properly rely on spatially and temporally distributed data in place of spatially aggregated and temporally distributed ones, such as those adopted in this work.

Conclusions
Exploring short-term climate change effects on vegetation through pheno-meteorological modeling and relative trends such as ET and PMs, snowpack time persistence to the ground derived from satellites EO data with different time scales provide useful information about biosphere-atmosphere exchanges of carbon, energy, and water at both regional and global scale. This information is crucial to calibrate adaptive strategies to better face climate change effects in alpine ecosystems. It is worth remembering that mountains are among the most sensitive ecosystems to climate change, especially rising temperature, where effects observed more quickly than in other terrestrial habitats [99][100][101][102][103][104][105][106][107][108][109][110][111][112][113][114] and this study could be considered another one concerning the short-term climate change effect on rangelands' and broad-lived forests' phenology in Aosta Valley. In particular, EO data from public and free archives of ready-to-use products such as MOD13Q1 and MOD16A2, proved to support well this type of analysis in spite of their reduced geometric resolution. Specifically, they proved to be useful in analyzing PMs and ET and SC time trends that, if compared with the air temperature trends, could support breeders, foresters and in general policymakers, to better address the management of alpine rangelands and broad-leaved woods. In this work, we showed that in Aosta Valley, phenological and evapotranspiration-related processes and snowpack melting time have been dramatically changed in the last two decades, with generally statistically significant trends. PMs, ET and snow cover time melting from EO data proved to significantly change their values in the last 20 years, with continuous progressive trends. Averagely in the area of study, rangelands and broadleaved forests, the length of season is getting longer, with a general advance in the SOS (start of the season) and a delay in the EOS (end of the season). With reference to ET, significant increasing trends were generally observed like other areas worldwide [114]. The water requirement from vegetation appeared to have averagely risen up by about Additionally, the snowpack timing persistence to the soil shows a general reduction trend through the years investigated. This pattern, together with an increase in air temperature anomalies, could be one of the main reasons that explain the PMs trends observed with a general boosting on the phenological season. The results encourage the adoption of remote sensing to monitor climate change effects on alpine vegetation, with a particular focus on the relationship between phenology and other abiotic factors, such as snowpack timing and temperature, supporting the idea that, presently, an effective technology transfer is now possible to dynamically support the management of mountain vegetation.