Exploring Climate Change Effects on Vegetation Phenology by MOD13Q1 Data: The Piemonte Region Case Study in the Period 2001–2019

: Rising temperature, rainfall, and wind regime changes, increasing of frequency and intensity of extreme events are only some of the effects of climate change affecting the agro-forestry sector. Earth Observation data from satellite missions (often available for free) can certainly support analysis of climate change effects on vegetation, making possible to improve land management in space and time. Within this context, the present work aims at investigating natural and agricultural vegetation, as mapped by Corine Land Cover (CLC) dataset, focusing on phenological metrics trends that can be possibly conditioned by the ongoing climate-change. The study area consists of the entire Piemonte region (NW-Italy). MOD13Q1-v6 dataset from TERRA MODIS mission was used to describe pluri-annual (2001–2019) phenological behavior of vegetation focusing on the following CLC classes: Non-irrigated arable land, Vineyards, Pastures, and Forests. After computing and mapping some phenological metrics as derivable from the interpretation of at-pixel level NDVI (Normalized Difference Vegetation Index) temporal proﬁle, we found that the most signiﬁcant one was the maximum annual NDVI (MaxNDVI). Consequently, its trend was analyzed at CLC class level for the whole Piemonte region. Natural and semi-natural vegetation classes (Pastures and Forests) were furtherly investigated testing signiﬁcance of the Percent Total Variation (TV%) of MaxNDVI in the period 2001–2019 for different altitude classes. Results proved that Non-irrigated arable land showed a not signiﬁcant trend of MaxNDVI; differently, vineyards and forests showed a signiﬁcant increasing one. Concerning TV%, it was found that it increases with altitude for the Forests CLC class, while it decreases with altitude for the pastures class.


Introduction
Humans have always interacted with environment; through their activities, they have modified and adapted natural environment by continuing to take advantage of resources over time. Starting from the second half of the 18th century, man has prevailed over environment and, during the last century, ecosystems structure has been changing much more rapidly than in the past [1]. According to the Global Environment Outlook [2], the planet is suffering the sixth great extinction that for the first time will be not due to natural global changes, but rather to human activity [3].
Rising temperatures, changes in rainfall and wind regime, and changes in frequency and intensity of extreme events are the main factors conditioning animal and plant species [4]. These events regulate environmental characteristics like the availability of nutrients, essential for the development of primary producers, ice cover and, in the sea, the intensity of convective movements, transparency, and water level. With these premises, it is mandatory to consider response times characterizing those processes influenced by climate change: from the short times for impacts on physiology (day-months) to the longer ones for area changes (year-decades), up to the typical scales of evolutionary processes (hundreds of year-millennia) [5,6].
Plants are very sensitive to climate and its variations: distribution of vegetative types defines the bioclimatic zones (or belts), whereas phenology is strictly dependent on seasonal meteorological trend [7,8]. Consequently, climate change could lead to changes in geographical distribution of species and, contemporarily, to timely shift appearance of plant pheno-phases, with probable consequences on crop productivity and landscape in general [9,10]. In particular, vegetation trend changes could be related to vegetation management, climate and climate change as deeply discussed in literature [11][12][13]. Phenological observations and agronomic calendars have been used by humans for thousands of years and long historical time series of such information can be found in many documents. There are, in fact, important historical archives: some records from Japan date back to the early 1800s [14]; some from Europe go back to the early 1700s [15]. More recently, over the last fifteen years, phenology has regained a new interest in the scientific community (included the European Topic Centre on Air and Climate Change (ETC-ACC) of the European Environment Agency (EEA)). Applications focusing on seasonal dynamics of vegetation [16][17][18] and phenological modelling use these data to predict capacity and productivity in agriculture [19] or as a source of robust proxies for medium and long-term effects of climatic variations.
Satellite data and predictive model simulations are the basis to understand the climate system [27]. Satellite remote sensing for Earth Observation can effectively support this process, making possible to get information about the earth's surface, subsoil, and atmosphere [28]. Since the first spatial observation of solar radiation and cloud reflection, remote sensing has gradually become a major research method in climate change studies [29]. The use of satellites permits the observation of states and processes of the atmosphere, earth, and oceans at different space-time scales. For example, it is one of the most efficient approaches for monitoring land cover and its changes over time at a variety of spatial scales [30,31]. Satellite data are often used in combination with climate models to simulate climate system dynamics and to improve climate forecasting [32]. The Global Climate Observing System (GCOS) has reported 26 of the 50 key climate variables (ECVs) as significantly dependent on satellite observations [33]. Remote sensing data are also largely used to develop prevention, mitigation, and adaptation measures to cope with the impact of climate change [34], to take out insurance policies based on vegetation trends [35][36][37][38], to support Common Agricultural Policy control in agriculture [39], to support wildlife diseases assessment [40], to manage the risk of falling trees and heat islands monitoring [41][42][43]. Studies on climate change, in fact, need continuous calibration and validation data and appropriate temporal and spatial sampling over a long period of time [44].

Study Aims
The aims of this work were to assess how low spatial resolution satellite data from TERRA MODIS mission can evaluate the effects of mid-term climate change on vegetated areas (that have not changed over time), with special concerns about annual phenological metrics as derivable from Normalized Difference Vegetation Index (NDVI) maps time series between 2001 and 2019 in the Piemonte Region.
A preliminary trend analysis was performed to test the presence of climate change in area of interest (AOI) according to regional meteorological stations. Consequently vegetated "natural/semi natural" (forests and pastures) and "managed" (agricultural) classes, derived by the Corine Land Cover (CLC), were, separately, considered. Moreover altitude effects on phenology changes were also taken into account to deepen investigation for forest and pasture classes.

Study Area
The study area corresponds to the entire Piemonte region (North-West Italy) ( Figure 1). It sizes about 25,400 km 2 . From a climatic point of view, North-Western Alps cause a gradual reduction of temperature while altitude increases making the local climate a typical temperate one with continental character. Average annual precipitation and temperature are 1050 mm and 11.9 • C, respectively. Thermal inversion phenomena, caused by cold air, can often affect the area, especially in the valleys. classes, derived by the Corine Land Cover (CLC), were, separately, considered. Moreover altitude effects on phenology changes were also taken into account to deepen investigation for forest and pasture classes.

Study Area
The study area corresponds to the entire Piemonte region (North-West Italy) ( Figure  1). It sizes about 25,400 km 2 . From a climatic point of view, North-Western Alps cause a gradual reduction of temperature while altitude increases making the local climate a typical temperate one with continental character. Average annual precipitation and temperature are 1050 mm and 11.9 °C, respectively. Thermal inversion phenomena, caused by cold air, can often affect the area, especially in the valleys.

Available Data
The following data were used:

Satellite Data
Monitoring of biomass along time requires that uniform and comparable metrics are used. Vegetation Indices (VIs) from multispectral data are often used to describe vegetation cover/status and, consequently, to assess phenology with high temporal resolution [45,46]. For this reason, data form the Moderate-Resolution Imaging Spectroradiometer (MODIS) sensor, operating on TERRA satellite mission since 2000 and having 1-day temporal resolution, was taken into account. The MOD13Q1-v6 NDVI collection product, obtained from NASA's Land Processes Distributed Active Archive Center (LPDAAC) [47], was specifically used. Data were downloaded from AppEEARS, provided as Tagged image File (TIF) georeferenced in the WGS84 geographic reference system. MOD13Q1-v6

Available Data
The following data were used:

Satellite Data
Monitoring of biomass along time requires that uniform and comparable metrics are used. Vegetation Indices (VIs) from multispectral data are often used to describe vegetation cover/status and, consequently, to assess phenology with high temporal resolution [45,46]. For this reason, data form the Moderate-Resolution Imaging Spectroradiometer (MODIS) sensor, operating on TERRA satellite mission since 2000 and having 1-day temporal resolution, was taken into account. The MOD13Q1-v6 NDVI collection product, obtained from NASA's Land Processes Distributed Active Archive Center (LPDAAC) [47], was specifically used. Data were downloaded from AppEEARS, provided as Tagged image File (TIF) georeferenced in the WGS84 geographic reference system. MOD13Q1-v6 product is a 16 days timely-spaced one, having a spatial resolution of 250 m. It is a composite product containing the "best" available local observations (at pixel level) within a 16 days period. Selection of the best pixel relies on criteria that minimize cloud cover effect and viewing angles and maximize NDVI local value (maximum in the reference period). MOD13Q1-V6 product contains many layers, included the NDVI and the pixel reliability (PR) ones. PR codes are reported in Table 1.  Table 2. CLC Level 3 is the most detailed one, where each patch presents more than 75% of the characteristics of a given class according to the CLC nomenclature rules [48]. All the available CLC maps were combined by intersection to recognize those patches whose meaning (class) remained the same along the entire explored period (2001-2019). Only agricultural and natural classes were considered (Table 3). Among these classes, only the Forests one includes sub-classes: deciduous forests (311), coniferous forests (312), and mixed deciduous and coniferous forests (313). Consequently, they were preventively merged into a single class (Forests) in order to define a single land cover class representative of natural vegetation conditions (no human activity is generally performed to manage tree growth). All data were finally projected into the WGS84/UTM 32N reference frame. Data processing was achieved by QGIS 3.2.0 and SAGA GIS 7.0. Statistical analysis was performed using Past 4.01.  [49]. It covers the entire regional territory and was generated from LiDAR (Light Detection and Ranging) point clouds. Native DTM, having a grid size of 5 m and a height accuracy of 0.66 m, was generalized down to a 25 m grid product, and made available for large size areas investigation. In this work, DTM was used to investigate dependence of natural and semi-natural vegetation (Pastures and Forests classes) from altitude; 4 altitudinal zones (hereinafter called AZs) were defined (Table 4) according to Mayr-Pavari bioclimatic zones [50]. In order to test climate induced trends, some ground-retrieved temperature data were used as reference. A total of 24 meteorological stations (MS) were selected from the Regional meteorological network (www.arpa.piemonte.it, accessed on: 16 January 2021) in order to explore the temperature trends in AOI. In particular daily average temperature was analyzed in the period 2000-2019 defining time series with a time step of 1 day (7671 observations for each MS). The MS selection was achieved by spatial intersection between the CLC and AZs classes previously defined. Specifically 3 MS temperature time series were collected for each CLC and AZs classes ( Figure 2).

Reference Classes and Patches Selection
Selection of unchanged CLC patches, i.e., maintaining the same class meaning along time, is mandatory while investigating climate trends along time (2001-2019) at class level. This step is basic to avoid misleading deductions that could significantly compromise time trend analysis. For example, NDVI trend of a forest area will certainly change if, in the meantime, it turns to urban. For this reason, a preliminary step was performed to select only those patches (belonging to the above-mentioned classes) whose meaning remained the same in all the available 4 CLC maps. Selection was achieved by polygons intersection by ordinary GIS procedures. Results are reported in Figure 3.
Since natural and semi-natural vegetation is assumed to be more susceptible by climate change [51] in respect of agricultural vegetation, a further level of investigation was carried out for natural classes solely (Forests and Pastures).

Reference Classes and Patches Selection
Selection of unchanged CLC patches, i.e., maintaining the same class meaning alon time, is mandatory while investigating climate trends along time (2001-2019) at clas level. This step is basic to avoid misleading deductions that could significantly compro mise time trend analysis. For example, NDVI trend of a forest area will certainly chang if, in the meantime, it turns to urban. For this reason, a preliminary step was performe to select only those patches (belonging to the above-mentioned classes) whose meanin remained the same in all the available 4 CLC maps. Selection was achieved by polygon intersection by ordinary GIS procedures. Results are reported in Figure 3. Since natural and semi-natural vegetation is assumed to be more susceptible by climate change [51] in respect of agricultural vegetation, a further level of investigation was carried out for natural classes solely (Forests and Pastures).
For these two classes authors investigated the role of altitude (at AZs level) on time For these two classes authors investigated the role of altitude (at AZs level) on time trends of the explored phenological metrics. CLC patches were consequently completed with the correspondent height information derived from DTM by GIS zonal statistics. Forest and Pasture classes were finally analyzed at AZs level determining a total of eight new sub-classes to analyze.
Two analyses were consequently carried out: one testing phenological metrics trend at class level with reference to classes of Table 3; one testing differences between phenological metrics trends of the same class placed at different AZ. For this second analysis Forest and Pasture classes, solely, (natural/semi-natural) were considered.

Temperature Trends Analysis
Starting from the MS data, the average time series was computed for each CLC and AZs classes. Resulting daily average temperature time series were analyzed by first order polynomial regression. In particular the gain values and their significance level were obtained in order to test the presence and quantify the strength of climate induced trends.
The basic assumption was that vegetation phenology is strongly affected by temperature [52,53]. Therefore an increase of temperature trend changes phenological response [54,55]. Nevertheless, human management can mitigate the climate change effects on vegetation [56].

Phenological Metrics (PM)
For each investigated class, the correspondent average NDVI temporal profile was computed from NTS (hereinafter called ANTS). Phenological response is strictly related to many factors as vegetation species, terrain aspect, local nutrients availability, and climate. Only the latter can be assumed homogenous over AOI while other factors are characterized by high variability necessarily affecting phenological metrics estimates. Averaging phenological profiles within the same land cover class mitigates this variability enhancing climate related effects [57,58].
The following phenological parameters were computed at year level with the aim of summarizing annual growing season (from leaf development to leaf senescence) of vegetated classes: SOS (Start Of Season), EOS (End Of Season), LOS (Length Of Season), MaxNDVI (annual maximum NDVI value), and DOY (Day Of the Year when MaxNDVI occurs). In particular SOS is define as the moment of the year having a consistent upward trend in leaf development. Otherwise the EOS is the moment of the year having a consistent downward trend in leaf senescence. While MaxNDVI is related to maximum biomass development and maximum level of canopy photosynthetic activity [60].
SOS and EOS were calculated by a first derivative approach where SOS and EOS are located at the maximum and minimum of the 1st derivative, respectively. These moments occur where ANTS is positively and negatively steeper [61,62]. LOS was obtained by differencing of EOS and SOS; it measures the vegetative season length in the considered year.
The previously mentioned metrics describe specific moments (in the time-domain) when vegetation expresses a peculiar event (e.g., leaf development or leaf senescence), but are not able to measure the strength of canopy vigor. Many authors reported that yearly maximum NDVI value is a good metric to compare vegetation behavior throughout the years [63,64]. Consequently, yearly MaxNDVI was calculated and the corresponding DOY recorded.

Class Effects on PM Trends
A feasibility analysis was conducted considering the 4 CLC classes of Table 3 to determine if and how phenological metrics showed significant trends, in the reference period, possibly relatable to climate change. All the above mentioned phenological metrics were consequently compared along the explored 19 years and fitted by a first order polynomial

Altitudinal Effects on PM Trends
A further assessment was done to refine deductions about natural (Forests) and seminatural (Pastures) vegetation. These 2 classes were split in 4 sub-classes each, according to the 4 defined AZs. Investigation concerned only the MaxNDVI metric, that the previous step indicated as the most significant one from the time trend point of view. Again, a first order polynomial linear regression was used to test strength and significance of trends. The approach is coherent with other similar works where NDVI-based vegetation changes and their responses to climate change over mountain system exploring different time range [29,66,67]. The adopted approach is the common slope test that is based on gain significance analysis [68,69] testing the hypothesis that two gain values are statistically different. It is worth to remind that gain value only defines the rate of increase/decrease of metric trends; therefore, to quantify the total variation that the considered PM suffered from in the reference period (2001-2019), the Percent Total Variation of MaxNDVI (TV%) was calculated according to Equation (1).
where G is the gain value of MaxNDVI linear regression model, O is the offset value and t 1 and t 19 are the first and the last observation respectively (2001, 2019). A final comparison to test dependency of TV% from AZs was performed by scatter plot analysis.

Reference Classes and Patches Selection
To map unchanged CLC classes during the period 2001-2019 a polygon layer intersection step was achieved with reference to the four CLC available maps. Statistics concerning focus classes (Non-irrigated arable land, Vineyards, Forests, and Pastures) are shown in Table 5. It can be noted that size of all the considered CLC classes showed a decreasing trend over time. In particular, pastures decreased of 1521 km 2 between 2006 and 2012, correspondent to about 70% of its value. Similar results were found by Rusu et al. in Romania for natural/semi natural areas related to deforestation and urbanization processes [70]. Future investigations could be relevant to discover the nature of this significant change in this particular area. Conversely, other classes remained almost unchanged along the entire period as also observed in other works [70]. As previously mentioned, forest and pasture classes were further divided into four sub-classes depending on AZs, whose dimension is reported in Table 6.  Table 6 shows that forests cover prevail at lower altitudes and decreases while altitude increases. More than 88% of forests class is, in fact, located between 0 and 1200 m a.s.l. Pasture class, differently, prevails at the higher altitudes (AZ4 > 2000 m a.s.l.). Table 7 shows gain values of linear regression trends according to CLC and AZs classes. All trends were significantly increased during the 2000-2019 period, suggesting the presence of climate change within AOI. Similar increasing trends were reported in literature [71,72], suggesting a key role of temperature in phenological development. Therefore authors used PM to assess the effects of such temperature increase (see Section 3.3). Table 8 shows gain values computed for all PM and CLC classes. In particular, not for all classes, the most significant PM was MaxNDVI. This metric seems to be more sensitive to vegetation changes along time series as proved by some works [61,73] that, similarly, adopted low geometrical resolution satellite imagery. Other PMs did not show significant trends for any CLC class. The same finding was also reported in literature [56,[74][75][76], in particular Tao [74] founds that MaxNDVI had a positive trend in croplands, forests and grasslands, more significant than LOS. This phenomenon could be related to two main factors: one related to the adopted product (MOD13Q1); one related to the PM itself. As far as the first one is concerned, it is worth to remind that MOD13Q1 product is a collection made of composited images, where pixel values refer to different dates, within a 16 days period. This, considerably, determines a significant uncertainty in SOS, EOS, and LOS estimation. As far as the second factor is concerned, PMs other than MaxNDVI are estimated looking for a single specific event along the analyzed NTS (e.g., the day of the end of season). Frequent harvests (or mowing) of crops, ordinarily performed in an agriculture context, certainly affect EOS/LOS estimates, especially when moderate resolution images are used, making them particularly uncertain. Among tested trends, only those related to Forests (natural) and Vineyards (agricultural with no water management) resulted significant. This suggests that human practices can mitigate effects of climate changes over vegetation; in fact, irrigated and fertilized agricultural classes did not show significant trend. Pastures is ordinarily considered a semi-natural system, but, unexpectedly, they did not show any significant trend; this is possibly due to their intensive management that is performed in Piemonte, especially in lowland areas. Similar results were obtained by Tang et al. [56] and Menzel et al. [56,71], specifically both studies showed that the effects of climate change on vegetation identified by PM changed in relation to the land cover. In particular, they noted that the agricultural class was poorly affected by climate change, while the forest one was the most susceptible to the climate effects.

Altitudinal Effects on PM Trends
MaxNDVI resulted the only PM showing a detectable and significant trend along time at the MOD13Q1 temporal and spatial scale. Consequently, a further analysis was performed to explore the dependence of its trends from altitude. Only natural and seminatural vegetation were considered at this point. Similarly, as previously done, first order polynomial trends were tested with reference to the eight classes of Table 6. Results are reported in Table 9 and Figure 4.  Table 10 reports p-values resulting from the common gain test performed for all the AZs of forest and pasture classes. In general, results indicate that no significant differences exist between gains of AZs trends for both forests and pastures. Nevertheless, according to Table 9, such trends show a positive value. These results suggest that climate change effects on MaxNDVI, i.e., vegetation vigor, is constant and latent in all the altitudinal classes. Nevertheless, this latent climatic-induced trend could vary in terms of relative strength if the vigor total variation (TV%, Equation (1)) is computed and compared with the initial MaxNDVI value.   TV% is reported against AZs in Figure 5 for both forests and pastures. Forests showed an increasing TV% gradient with AZs. Specifically, lower TV% (about 1%) correspond to forests placed at lower altitudes, suggesting that climate change determine stronger effects on forests located at higher altitudes (about 4%). Results seem to be consistent with a recent studies [74,[77][78][79] that found that vegetation greening and climate warming effects trends were stronger at higher altitudes, thus supporting the idea that an altitudinal gradient effect operates in determining forests vigor level and changes.
TV % is reported against AZs in Figure 5 for both forests and pastures. Forests showed an increasing TV % gradient with AZs. Specifically, lower TV % (about 1 %) correspond to forests placed at lower altitudes, suggesting that climate change determine stronger effects on forests located at higher altitudes (about 4 %). Results seem to be consistent with a recent studies [74,[77][78][79] that found that vegetation greening and climate warming effects trends were stronger at higher altitudes, thus supporting the idea that an altitudinal gradient effect operates in determining forests vigor level and changes.  Conversely, pastures showed a negative altitudinal TV% gradient. Pastures located at lower altitudes were more affected by climate changes than high altitude pastures. This probably relies on the fact that pastures, in lowlands, are commonly more intensively managed (irrigations and fertilizations) than in mountain. In this, human management could increase vigor and productivity of low altitude grass systems while grass systems in alpine zone (i.e., high AZs) are generally characterized by low biomass and often alternate with bare soil [80,81]; this factor certainly affects ANTS determining a minor sensitivity of phenology to climate changes. The same phenomenon was reported in Li [82] where a significant NDVI increase was primarily associated with forests, while a weak NDVI increase was observed for pastures at altitudes > 2000 m. Results are also in agreement with [83] findings that proved that main vegetation climate-related changes occur at lower altitudes for alpine pastures.

Conclusions
In this work, the increasing trend of temperature probably related to climate change in AOI was proved analyzing regional meteorological stations data. Therefore effects on vegetation phenology was expected. To prove such phenomenon a trend analysis of climate change effects on Piemonte Region vegetation was performed using MOD13Q1 imagery in the period 2001-2019. In particular, agricultural and natural land cover classes (as reported by CLC maps) were explored by analyzing correspondent ANTS. Some PMs were computed with reference to the yearly class ANTS, but only MaxNDVI showed significant trends possibly related to climate change for some classes at regional level. In particular, agriculture classes did not show any significant MaxNDVI trend (G = 11.285, p > 0.05); differently, forests and vineyards presented a significant positive trend (G = 10.500, p < 0.05 and G = 34.333, p < 0.01 respectively). A further analysis was conducted to assess if altitude could condition trends of natural and semi-natural classes. The common gain test proved that no significant difference exists between trends of forests classes placed at different altitudes. Similarly, for pastures. Nevertheless, while relating the absolute difference of MaxNDVI totally occurred along the entire period (2001-2019) with its initial value at 2001 (by TV%), it was observed that high altitude forests presented a higher value of TV% (4%) than low altitude ones (1%). Conversely, pastures showed higher TV% at lower (7.7%) than at higher altitudes (3.3%). In conclusion, MOD13Q1 data proved to be an effective product to monitor vegetation over time and to derived some PMs. Unfortunately, not all of them, at the temporal and spatial scale of MOD13Q1 and at the regional level of class aggregation proved to be sensitive to the ongoing climate change. Only MaxNDVI showed significant trends, at regional level, for forests and vineyards making evident that, absolutely speaking, it operated similarly at different altitude zones within the same class. Nevertheless, we showed that same absolute difference in MaxNDVI meant a different effect (measured by TV%) at different altitudes, depending on its initial value.