NDVI or PPI: A (Quick) Comparison for Vegetation Dynamics Monitoring in Mountainous Area

Abstract

Cold ecosystems are experiencing a warming rate that is twice as fast as the global average and are particularly vulnerable to the consequences of climate change. In mountain ecosystems, it is particularly important to monitor vegetation to understand ecosystem dynamics, biodiversity conservation, and the resilience of these fragile ecosystems to global change. Hence, we used satellite data acquired by Sentinel-2 to perform a comparative assessment of the Normalized Difference Vegetation Index (NDVI) and the Plant Phenology Index (PPI) in mountainous regions (canton of Valais-Switzerland in the European Alps) for monitoring vegetation dynamics of four types: deciduous trees, coniferous trees, grasslands, and shrublands. Results indicate that the NDVI is particularly noisy in the seasonal cycle at the beginning/end of the snow season and for coniferous trees, which is consistent with its known snow sensitivity issue and difficulties in retrieving signal variation in dense and evergreen vegetation. The PPI seems to deal with these problems but tends to overestimate peak values, which could be attributed to its logarithmic formula and derived high sensitivity to variations in near-infrared (NIR) and red reflectance during the peak growing season. Concerning seasonal parameters retrieval, we find close concordance in the results for the start of season (SOS) and end of season (EOS) between indices, except for coniferous trees. Peak of season (POS) results exhibit important differences between the indices. Our findings suggest that PPI is a robust remote sensed index for vegetation monitoring in seasonal snow-covered and complex mountain environments.

1. Introduction

Cold ecosystems (i.e., arctic and alpine regions) are affected by a warming rate at least twice higher than the global average and, therefore, are particularly sensitive to climate change impacts [1,2,3]. Temperature [4,5,6,7] and snow cover duration [7,8,9,10] are known as major limiting factors for vegetation growth in these environments, and both are influenced by climate change.

Monitoring vegetation productivity and phenology in mountainous areas is essential for understanding ecosystem dynamics, biodiversity conservation, and the resilience of these fragile ecosystems to global change [11,12,13,14]. Vegetation plays a key role in mountain ecosystems by providing important services such as habitat for a wide range of plant and animal species, grazing resources for livestock, regulating water and nutrient cycles (e.g., quality and quantity), stabilizing soil against erosion or offering recreational opportunities for tourism and outdoor activities [15,16,17,18]. However, these ecosystems are particularly vulnerable to climate or land-use/cover (LUC) change [8,19]. Monitoring changes in vegetation cover, species composition, and ecosystem structure over time can help identify emerging threats and guide effective land management strategies, conservation measures, and restoration efforts to maintain the ecological integrity and functioning of mountain landscapes [7].

In the Alps, the vegetation seasonal phenological cycle occurs during the snow-free period [20]. Vegetation phenology refers to the study of recurring events in the life cycle of plants (i.e., budding, flowering, senescence) reflecting the dynamics of terrestrial ecosystems [21]. Vegetation phenology is particularly sensitive to climate change as well as land-cover and land-use change [21], inducing various consequences such as phenological shifts (e.g., change in start, duration, end of season) or change in vegetation productivity [6,8,21,22], strongly affecting mountainous ecosystems [12,23]. As mentioned by [24], the length of the growing period is one of the main determining factors of productivity. In addition, a global vegetation greening trend has been identified since the 1980s [25,26]. In the cold areas, the greening is consistent with the recent temperature increase [26]. Although this trend is well-documented in the Arctic [4,27], there are only a few studies in the European Alps [8,11,28]. However, there is an incomplete understanding of the mechanistic links between climate change and these changes in vegetation dynamics [27,29].

Recent technological advancements, in particular remote sensing, offer new possibilities for LUC change [30,31,32,33] and vegetation dynamics monitoring, particularly in wide and complex ecosystems [21,34]. Remote sensing and derived vegetation indices offer the opportunity to observe vegetation dynamics over large areas and long periods [5,34,35]. Depending on their spatial, temporal, and spectral resolutions, different sensors may provide various opportunities for vegetation monitoring [35]. For global, continental, and national mapping, the MODerate resolution Imaging Spectroradiometer (MODIS) is extensively used [8,36,37] due to its 1 to 2 days temporal resolution and 250 to 1000 m spatial resolution (depending on the spectral band). Landsat satellites are widely used on a regional to national scale with a spatial resolution of 30 m and a temporal resolution of 16 days [4,11,38]. More recently, Sentinel-2 satellite constellations and their multispectral products offer new observation possibilities for continental to local mapping with a spatial resolution of 10 m and a revisit time of 5 days [39], providing new capabilities for more detailed ecosystems monitoring, such as mountain vegetation dynamics or phenological metrics retrieval at intra- and inter-annual time steps [34].

Among the different satellite-derived vegetation indices, such as the Enhanced Vegetation Index (EVI), the Green Normalized Difference Vegetation Index (GNDVI) the kernel Normalized Difference Vegetation Index (kNDVI), the Chlorophyll Vegetation Index (CVI,) or the Plant Phenology Index (PPI), the Normalized Difference Vegetation Index (NDVI) is the most widely used index for assessing vegetation dynamics, specifically in mountain environments [8,11,35,39,40,41]. NDVI measures the difference between near-infrared (NIR) and red light reflectance, providing information on vegetation greenness and photosynthetic activity. It is widely used to monitor vegetation density, health, and productivity over large areas. However, it suffers from various limitations for vegetation state retrieval, particularly in seasonal snow-covered environments and coniferous forests [36]. The Plant Phenology Index (PPI) has been recently developed to overcome these limitations and improve plant phenology monitoring, notably in Arctic areas, and has shown promising performances [40,42]. Unlike the NDVI, the PPI is physically based and derived from Beer’s law, adapted for canopy reflectance. It has an almost linear relationship with the Leaf Area Index (LAI), making it a powerful proxy for phenological stage identification through remote sensing [40]. Recent studies have further explored the potential of PPI beyond Arctic regions, offering new insights into vegetation monitoring in arid ecosystems [43] and across Europe [39]. To our knowledge, the PPI performance has not been assessed specifically for complex mountainous environments.

As in Arctic environments, mountain areas, and more specifically, the European Alps, face a seasonality marked by the presence of snow and low temperatures. Therefore, the PPI could provide new insights for coping with the limitations of NDVI and enhance the monitoring of vegetation in complex environments such as mountain areas. Considering these factors, the objectives of this paper are to investigate (1) how PPI compares to NDVI in retrieving phenological metrics in areas characterized by seasonal snow and coniferous forests, such as mountainous regions, and (2) whether PPI can address some of the known limitations of NDVI for monitoring vegetation dynamics in such environments. These aspects will be examined through an initial comparative assessment of phenological metrics derived from PPI and NDVI for monitoring vegetation dynamics across four vegetation types—deciduous trees, coniferous trees, grasslands, and shrublands—in European alpine ecosystems. The findings of this study may provide new perspectives on monitoring vegetation dynamics in complex environments, potentially linking long-term regional vegetation trends to annual phenological and growth responses to rising temperatures and climate stress in diverse topoclimatic conditions. Ultimately, this could aid in forecasting changes in ecosystem productivity and biodiversity.

2. Materials and Methods

2.1. Study Area

The study area is delimited through the borders of the canton of Valais in Switzerland (Figure 1). The lowest point is at 372 m at Geneva Lake, and the highest is at 4634 m at the Dufourspitze. The surface of the area of interest is about 5224.35 km². The Valais encompasses several valleys, with the Rhône Valley being the central one. This mountainous region is situated in the northwestern European Alps, with several peaks rising above 4000 m. This mountain range shelters the inner Rhône Valley from moist oceanic air masses, resulting in a steep precipitation gradient from the west (around 1000 mm/year) to the east (around 600 mm/year) [43,44,45]. The eastern Rhône Valley is considered one of the driest and most continental inner alpine valleys and hosts a rich xerothermic flora [44]. Cloud cover is generally low, while solar radiation is high [45]. As in other mountainous regions, the Valais is therefore characterized by a high diversity of vegetation types and ecosystems due to its diverse geographical features [7].

2.2. Methodology and Implementation

The steps of the implemented methodology are shown in Figure 2. We used the R environment (https://www.R-project.org/ accessed on 3 August 2023) for image processing and statistical analysis, TIMESAT 3.3 for time-series processing and phenological metrics retrieval [46,47], and Google Earth Engine (GEE) [48] and ArcGIS pro 3.1 software for vegetation cover mapping.

2.2.1. Satellite Data Access and Processing

Satellite data have been gathered for the vegetation cover classification and the comparison of PPI and NDVI over the Valais territory.

For the vegetation cover classification, we used the Forest Type (https://land.copernicus.eu/en/products/high-resolution-layer-forest-type/forest-type-2018 accessed on 17 August 2023) and Grassland (https://land.copernicus.eu/en/products/high-resolution-layer-grassland/grassland-2018 accessed on 17 August 2023) high-resolution layers of the Copernicus Land Monitoring Service (CLMS). These data have a resolution of 10 m, characterizing broadleaf forest, coniferous forest, and grasslands. For mapping shrublands, we followed the methodology proposed by [49], using Sentinel-2 data in the Google Earth Engine to compute first the median of the Normalized Anthocyanin Reflectance Index (NARI) and the Normalized Chlorophyll Reflectance Index (NCRI) for each year between 2017 and 2021 from September 1st to November 1st and then computed the median of the five years. The detailed methodology is explained in the next section.

All PPI and NDVI Sentinel-2 scenes for the years 2018 to 2022 were downloaded from the WEkEO platform (https://www.wekeo.eu/ accessed on 10 August 2023), corresponding to 4 tiles for covering the Valais region: T32TLS, T32TMS, T32TMR, T32TLR. There are 73 observations per year at 5-day intervals, beginning on 05-01-2018, resulting in 1460 images per index, 2920 in total. In addition, for each observation date, a quality layer (i.e., QFLAG2) was downloaded. It corresponds to a bitwise encoded status map, which indicates if the pixel is water, land, cloud, snow, or shadow [34].

2.2.2. Vegetation Cover

Regarding the vegetation cover classification, three vegetation classes have been extracted from the High-Resolution Layers of the CLMS: deciduous (broadleaved) and coniferous trees and grasslands for Valais. The limits of the canton have been obtained from the swissBOUNDARIES3D product (https://www.swisstopo.admin.ch/en/landscape-model-swissboundaries3d accessed on 03 August 2023).

To retrieve the shrublands layer, we followed the methodology proposed by [49,50]. All available Sentinel-2 scenes over Valais for the years 2017 to 2021 from September 1st to November 1st have been extracted in GEE. Then, the Sentinel-2 cloud probability product has been applied on every scene with a threshold of 0.65 to remove the clouds and clouds shadows. Third, the NARI and NCRI were computed on all scenes, and a median composite per year was generated. Finally, we computed a median of the 5 years.

The NARI is sensitive to the plant canopy anthocyanin content. Using this index, ref. [50] developed a methodology to improve the mapping of mountain shrublands. As they noted, shrublands are dominated in the European Alps by Ericaceae (i.e., Vaccinium spp. and Rhododendron ferrugineum). These species have the particularity to accumulate red anthocyanin pigments in the late autumn. This characteristic offers an opportunity to differentiate Ericaceae shrublands from other vegetation types. Indeed, as demonstrated by [50], their presence is often underestimated and confounded with grasslands. Therefore, Sentinel-2 data could be useful to enhance the discrimination of shrublands using the NARI index that is computed as follows [50]:

$$\mathrm{N}\mathrm{A}\mathrm{R}\mathrm{I}={\displaystyle \frac{\frac{1}{\mathsf{\rho } \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}-\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}}{\frac{1}{\mathsf{\rho } \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}}}$$

(1)

where ρ Green is the reflectance in the Green band, ρ Red-edge is the reflectance of the Red-edge band.

The NCRI is a normalized adjustment of the canopy chlorophyll content proposed by [49]. It allows us, in complement to NARI, to discriminate forest from Ericaceae shrublands and grasslands. NCRI is calculated using Sentinel-2 band 5 (Red-edge) and band 7 (Red-edge 2) [49]:

$$\mathrm{N}\mathrm{C}\mathrm{R}\mathrm{I}={\displaystyle \frac{\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}-\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}2}}{\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}+\frac{1}{\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}2}}}$$

(2)

where ρ is the reflectance in the respective Red-edge or Red-edge 2 band.

Following the results of [49], we applied a threshold of >0.325 for the NARI and a threshold of <0.42 for the NCRI to discriminate among the Ericaceae shrublands.

The resulting shrublands layer has a resolution of 20 × 20 m. This raster layer has been resampled and reprojected using the nearest neighbor method to correspond to the 10 × 10 m resolution of the Copernicus High Resolution layers. It is important to note that this methodology does not allow for the identification of non-Ericaceae evergreen shrublands like Juniperus and, consequently, results in an incomplete mapping of shrublands [49,50].

The final vegetation cover map was obtained by attributing the corresponding classes to each pixel of Valais in a single raster layer.

2.2.3. Vegetation Indices

For vegetation indices images pre-processing, we first mosaicked the four tiles and extracted the value for our region of interest for PPI, NDVI, and Quality flag (QFLAG2). Second, following the recommendation of [34], we applied a medium filter by masking the overlying PPI and NDVI pixel for the QFLAG2 values ranging from 4 to 2048, corresponding to snow, cloud, and shadow filtering with the addition of surrounding pixels.

NDVI

As mentioned before, the NDVI is one of the most widely used remote sensing indices for vegetation monitoring. As noted by [40], the popularity of this index is due to its relative robustness against noise and sun-sensor geometry variations and the availability of long-term time-series at a global scale. The NDVI is calculated following:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{\mathsf{\rho } \mathrm{N}\mathrm{I}\mathrm{R}-\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}}{\mathsf{\rho } \mathrm{N}\mathrm{I}\mathrm{R}+\mathsf{\rho } \mathrm{R}\mathrm{e}\mathrm{d}}}$$

(3)

where ρ NIR is the reflectance in the near-infrared band and ρ Red is the reflectance of the red band. The result is used as a proxy of the land surface greenness [8,41]: a value of −1 indicates a water surface; a value between 0 and 0.2 corresponds to almost non-vegetated areas; a value close to 1 represents a dense and green vegetation cover. A healthy vegetation cover tends to have a high absorption of photosynthetically active radiation, captured in the red band, and a low absorption of low infrared radiation, which could induce damage to plants due to overheating [51]. This index, however, suffers from two major limitations [40]: (1) its sensitivity to soil background (e.g., snow) and (2) its saturation at high vegetation density. In particular, this index shows difficulties in adequately identifying small amplitude variations, like evergreen coniferous forests [39].

PPI

Ref. [40] have developed a physically-based spectral index to characterize the phenological dynamics of vegetation: the Plant Phenology Index (PPI). This index has an almost linear relationship with the green Leaf Area Index (LAI) and is based on the Near-Infrared (NIR) and Red reflectance values. It follows Beer’s law, modified for canopy reflectance. PPI has the same unit as LAI (m²∙m⁻²) and is formulated as [39,40,42,52]:

$$\mathrm{P}\mathrm{P}\mathrm{I}=-K\ast \mathrm{l}\mathrm{n}\left({\displaystyle \frac{\mathrm{D}\mathrm{V}\mathrm{I}max-\mathrm{D}\mathrm{V}\mathrm{I}}{\mathrm{D}\mathrm{V}\mathrm{I}max-\mathrm{D}\mathrm{V}\mathrm{I}s}}\right)$$

(4)

where DVI (Difference Vegetation Index) is the difference between NIR and Red reflectance; DVI_max is the maximum canopy DVI of a specific site; DVI_s is the soil DVI. K is a gain factor given by:

$$\mathrm{K}={\displaystyle \frac{1}{4\ast (dc+0.5\ast (1-dc)/\mathrm{c}\mathrm{o}\mathrm{s}(\mathsf{\theta }s)}}\mathrm{*}{\displaystyle \frac{1+\mathrm{D}\mathrm{V}\mathrm{I}max}{1-\mathrm{D}\mathrm{V}\mathrm{I}max}}$$

(5)

where dc is an instantaneous diffuse fraction of solar radiation when the sun is at zenith angle θs (obtained from the corresponding scene metadata), calculated as:

$$\mathrm{d}\mathrm{c}=0.0336+0.0477/\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathsf{\theta }s\right)$$

(6)

For further information about the PPI formulation, we refer to the reference paper of [40].

2.2.4. PPI and NDVI Time-Series for Different Vegetation Classes

As noted by [53], time-series of vegetation indices derived from satellite spectral data are instrumental in monitoring seasonal vegetation development both within a single season and across multiple years. To facilitate this analysis, they developed the TIMESAT software package [46,53,54]. TIMESAT is particularly effective for analyzing satellite sensor data time-series, enabling detailed exploration and retrieval of seasonality parameters such as the start, peak, and end of the growing season.

To obtain the PPI and NDVI values for each vegetation class and date, we first resampled the vegetation cover map using the nearest neighbor method to match the 10 m resolution and extent of PPI and NDVI images. Subsequently, we employed the zonal function of the terra R package [55] to extract the mean PPI and NDVI values for each vegetation type across all observations.

To further analyze vegetation indices (VIs), we constructed time-series and fitted double logistic functions using the TIMESAT 3.3 software [46,47]. Following the methodologies outlined by [53] and other studies employing TIMESAT [37,56], we first removed spikes and outliers by applying a median filter with a window size of 1.5. Next, we applied consistent data ranges for each VI: −1 to 1 for NDVI and 0 to 3 for PPI. Given that noise in VI data from remote sensing is typically negatively biased, we adapted the upper envelope with a strength of 3 to reduce the influence of low data values. Finally, we set the season parameter to 1 to fit one season per year. As [53] noted, the success of fitting functions to time-series is as much an art as a science, relying heavily on visual examination and the specific nature of noise and disturbances in the data.

2.2.5. Seasonality Parameters

We also used the TIMESAT software to extract seasonal parameters such as the start of season (SOS), peak of season (POS), and end of season (EOS). This was achieved by employing the method based on seasonal amplitude as described by [53]. Accordingly, SOS and EOS were determined when the fitted curve reached a specific fraction of the difference between the base level and the POS. Following the PPI calibration by [39], we applied thresholds of 0.4 for NDVI and 0.25 for PPI to determine SOS and thresholds of 0.5 for NDVI and 0.15 for PPI to determine EOS.

3. Results

3.1. Vegetation Cover Map

The resulting vegetation cover map (Figure 3) and the derived surface and altitude metrics (

Table 1. Vegetation classes with corresponding area, proportion of Valais territory, and altitude mean.

Land Cover Classes	Area [km2]	Proportion [%]	Mean Altitude [m]
Deciduous trees	484	9.3	1387
Coniferous trees	722	13.8	1558
Grasslands	1176	22.5	2093
Shrublands	133	2.5	2103
Other	2713	51.9	2444

) show specific vegetation patterns for the Valais. Deciduous trees are concentrated in valleys at lower altitudes and correspond to a proportion of 9.3% (i.e., 484 km²) of the entire Valais area. Coniferous trees are mainly located at higher altitudes and cover 13.8% of the canton territory, which corresponds to 722 km². Grasslands are the most widespread vegetation class with an 1176 km² surface area, which is equivalent to 22.5% of Valais. They are mainly located above the tree line, although they are also present in valleys. Their mean altitude is above 2000 m. Shrublands are the vegetation class with the higher altitude mean (i.e., 2103 m). They also have the smallest surface area with a proportion of 2.5% (i.e., 133 km²) of the entire Valais. Altogether, these four vegetation classes represent 49.1% of the canton, the rest being mainly bare rocks, water, and built-up areas.

3.2. Vegetation Indices Time-Series

The raw NDVI time-series (Figure 4) are particularly noisy during the winter and spring, especially for coniferous trees. For deciduous trees, the metrics derived from the double logistic function indicate a mean NDVI of 0.54 over the entire period, with minimum and maximum values of 0.22 and 0.82, respectively. Coniferous trees exhibit a mean NDVI of 0.59, with corresponding minimum and maximum values of 0.23 and 0.80. Grasslands show a lower mean NDVI of 0.34, with a range from 0.06 to 0.56. Shrublands have a slightly higher mean NDVI of 0.40, with minimum and maximum values of 0.04 and 0.77. All vegetation classes exhibit a marked seasonal cycle.

For the PPI results (Figure 5), the raw data also display noise, albeit to a lesser extent, particularly in the winter and spring. The metrics derived from double logistic functions show mean PPI values of 0.52 for deciduous trees, 0.35 for coniferous trees, 0.43 for grasslands, and 0.41 for shrublands. The paired minimum and maximum values (min/max) are 0.01/1.50 for deciduous trees, 0.06/0.96 for coniferous trees, 0.01/1.13 for grasslands, and 0.01/1.54 for shrublands.

3.3. Time-Series: By Vegetation Class

In comparing the time-series of vegetation indices (VIs), we observe that the PPI exhibits a higher amplitude than the NDVI across all vegetation classes (Figure 6). This difference is less pronounced for coniferous trees. Additionally, PPI peaks are more distinct for deciduous trees, grasslands, and shrublands. The annual shape of the indices also differs, with PPI showing a tighter and more elongated pattern. For coniferous trees, PPI and NDVI display different scaled values.

A Shapiro test [57] revealed that the time-series data are not normally distributed. Consequently, we performed a Spearman correlation test [58] for each pair of vegetation classes [28,38]. The Spearman correlation coefficients (rho) are as follows: 0.91 for deciduous trees, 0.87 for coniferous trees, 0.93 for grasslands, and 0.96 for shrublands, each with a p-value less than 0.05.

3.4. Seasonality Parameters Retrieval

The seasonality parameters indicate close concordance between the two indices (PPI and NDVI) in determining the SOS and EOS for deciduous trees, grasslands, and shrublands (

Table 2. Seasonality parameters per vegetation type in day of the year (DOY) for PPI and NDVI and in day for differences between PPI and NDVI.

LC Classes		PPI				NDVI				Difference
Deciduous trees	Year	SOS	EOS	Length	POS	SOS	EOS	Length	POS	SOS	EOS	Length	POS
	2018	113	303	190	183	115	310	195	214	−2	−7	−5	−31
	2019	125	309	184	198	118	297	179	218	7	12	5	−20
	2020	110	295	185	184	109	291	182	197	1	4	3	−13
	2021	129	306	177	196	121	308	187	232	8	−2	−10	−36
	2022	120	301	181	181	127	315	188	223	−7	−14	−7	−42
	Mean	119.4	302.8	183.4	188.4	118	304.2	186.2	216.8	1.4	−1.4	−2.8	−28.4
Coniferous trees	Year	SOS	EOS	Length	POS	SOS	EOS	Length	POS	SOS	EOS	Length	POS
	2018	121	291	170	188	128	318	190	233	−7	−27	−20	−45
	2019	129	307	178	202	92	292	200	221	37	15	−22	−19
	2020	107	300	193	189	88	317	229	205	19	−17	−36	−16
	2021	124	299	175	198	97	311	214	232	27	−12	−39	−34
	2022	126	299	173	186	38	319	281	216	88	−20	−108	−30
	Mean	121.4	299.2	177.8	192.6	88.6	311.4	222.8	221.4	32.8	−12.2	−45	−28.8
Grasslands	Year	SOS	EOS	Length	POS	SOS	EOS	Length	POS	SOS	EOS	Length	POS
	2018	103	303	200	190	111	309	198	217	−8	−6	2	−27
	2019	104	303	199	203	101	297	196	218	3	6	3	−15
	2020	99	285	186	192	91	301	210	200	8	−16	−24	−8
	2021	96	301	205	203	88	307	219	216	8	−6	−14	−13
	2022	108	291	183	184	112	318	206	217	−4	−27	−23	−33
	Mean	102	296.6	194.6	194.4	100.6	306.4	205.8	213.6	1.4	−9.8	−11.2	−19.2
Shrublands	Year	SOS	EOS	Length	POS	SOS	EOS	Length	POS	SOS	EOS	Length	POS
	2018	147	294	147	203	144	301	157	226	3	−7	−10	−23
	2019	158	295	137	212	157	299	142	228	1	−4	−5	−16
	2020	138	288	150	204	120	309	189	214	18	−21	−39	−10
	2021	159	291	132	217	158	296	138	230	1	−5	−6	−13
	2022	145	285	140	195	147	307	160	223	−2	−22	−20	−28
	Mean	149.4	290.6	141.2	206.2	145.2	302.4	157.2	224.2	4.2	−11.8	−16	−18

). The mean differences (PPI−NDVI) for the SOS are 1.4 days for deciduous trees and grasslands and 4.2 days for shrublands. For the EOS, the mean differences are −1.4 days for deciduous trees, −9.8 days for grasslands, and −11.8 days for shrublands. However, for coniferous trees, the differences are more significant, with a mean difference of 32.8 days for the SOS and −12.2 days for the EOS. Regarding the length of the season, the divergences are minimal for deciduous trees (−2.8 days), followed by grasslands (−11.2 days), shrublands (−16 days), and coniferous trees (−45 days). For the POS retrieval, significant differences are observed across all vegetation classes, with PPI detecting the peak earlier than NDVI.

4. Discussion

The results of this study highlight the potential of the PPI in monitoring vegetation phenology in mountainous regions. The comparison of PPI and NDVI time-series depicts a less noisy signal for PPI, particularly during the transition seasons when snow could be present and for coniferous trees. For the retrieval of seasonality parameters, we found close concordance between the two indices for the SOS, EOS, and length of season for deciduous trees, grasslands, and shrublands. However, this is not the case for POS and for all seasonality parameters for coniferous trees. In the following section, we will discuss these results relying on existing literature.

4.1. Time-Series of the PPI and NDVI in Mountainous Areas for Different Vegetation Types

Our results suggest that the PPI is a robust remotely sensed proxy for time-series retrieval and phenological monitoring, which could complement or substitute for NDVI in mountainous areas. Although NDVI is widely recognized for its ease of processing, it has notable limitations in snow-covered environments and dense vegetation canopies such as evergreen forests [36,39,40]. These limitations are evident in the pronounced noise during winter and at the beginning and end of the growing season (Figure 4). The length of the growing season is a critical determinant of vegetation productivity [24], making precise retrieval of the SOS and EOS crucial. NDVI’s sensitivity to snow is also characterized by high variability at the start and end of the snow season [40]. Additionally, NDVI tends to saturate coniferous trees at levels similar to deciduous trees (approximately 0.8, Figure 4 and Figure 6). [39] reported that NDVI struggles to accurately retrieve Gross Primary Production (GPP) for coniferous trees, resulting in a negative correlation for the sites studied.

The PPI appears to mitigate these issues (Figure 5). However, consistent with previous studies [39,40], the PPI tends to accentuate the signal during the peak growing season (Figure 5 and Figure 6). This is due to its logarithmic formula, which increases sensitivity to variations in NIR and red reflectance [39].

4.2. Seasonality Parameters Retrieval

By using distinct thresholds for seasonality parameter retrieval for both the PPI and NDVI, we observed close concordance between the two indices for the SOS, EOS, and length of the season retrieval for deciduous trees, grasslands, and shrublands (Table 2). However, this concordance does not extend to the POS and all parameters for coniferous trees.

For SOS and EOS retrieval, we adopted a methodology based on amplitude thresholds identified by [39], which demonstrated the best alignment between PPI, NDVI, and ground-observed phenological stages. This approach yielded similar results between the two indices (except for conifers). However, [42] previously illustrated in a study aimed at distinguishing remotely measured plant phenology and snow seasonality that NDVI-derived land surface phenology correlates better with snow seasonality than with actual plant phenology, in contrast to PPI, which aligns well with ground phenology observation and Gross Primary Production (GPP) dynamics. This alignment is attributed to the linearity between PPI and Leaf Area Index (LAI), as proposed by [39]. They argue that the performance of Vegetation Indices (VIs) is not solely determined by vegetation change identification but also by their robustness against background noise. Overall, the use of PPI demonstrates superior performance for SOS and EOS retrieval in areas with seasonal snow cover and should, therefore, be prioritized [39,42,56].

Regarding the discrepancies in POS identification, the tendency of PPI to overestimate seasonal peaks likely exerts an influence. [39] suggest the possibility of applying an outlier filter to address this issue. For coniferous trees, the high level of noise in the signal and the plateau observed in NDVI (Figure 4 and Figure 6) are likely attributable to its known difficulties in retrieving the signal for this vegetation class, consequently influencing all seasonality parameter retrieval [56]. Therefore, PPI appears to be more reliable for dense coniferous forests as well.

4.3. Limitations

Enhancing the accuracy (i.e., minimizing the effects of atmospheric perturbation) of this study could be achieved by utilizing Level-2A Sentinel-2 images (i.e., surface reflectance) instead of Level-1C (top of atmosphere) imagery. To implement this, we recommend employing the Sen2cor program [59] to correct atmospheric effects in Level-1C images. The preprocessing of the substantial volume of images utilized in this study (n = 2920) posed challenges in terms of computational time and resources. For future investigations, such as pixel-scaled spatial analysis, considering high-performance computing is advisable. Recent advancements in data cubes for efficiently processing large volumes of satellite analysis-ready data offer a promising avenue for reducing preprocessing time, thereby enabling more comprehensive studies [60,61,62,63].

It is pertinent to note that our vegetation mapping method does not facilitate the identification of evergreen non-Ericaceae shrublands, leading to an incomplete characterization of this vegetation class and its dynamics. To address this limitation, a time-first approach for land cover mapping could offer significant advantages [32]. By retrieving dense time-series data of Sentinel-2 imagery and calculating various indices to monitor the seasonal spectral signals of different types of vegetation, it becomes possible to discriminate them [64].

Furthermore, while our study primarily focuses on the behavior of four main vegetation classes across the entire Valais region, there are avenues for further exploration. It would be beneficial to (i) incorporate more detailed vegetation types and (ii) evaluate spatial trends at pixel resolution to gain deeper insights into the vegetation dynamics.

It would also be highly valuable to analyze more in-depth relationships with ground data, such as eddy covariance GPP derived from FluxTower [65,66], for several locations across the Alps. [39] demonstrated that PPI is closely related to eddy covariance GPP data throughout Europe. Unfortunately, there are no data specifically available for the Valais. The nearest FluxTower is in Torgnon, in the Val d’Aoste (IT). This underscores the need for more high-quality, spatially distributed ground data in the Alps.

4.4. Contributions and Perspectives

The primary contribution of this study is to explore the potential of PPI in alpine regions, characterized notably by seasonal snow cover and the abundance of coniferous forests. The results are promising, indicating that PPI offers a reliable opportunity for remote sensing-based phenological and productivity monitoring in these challenging environments. Furthermore, our study underscores the significance of the recently launched Sentinel-2 satellite constellation, which provides enhanced possibilities for more accurate studies with its 10 m resolution and 5-day revisit time, particularly in complex terrain such as mountainous areas. To advance this research field, we recommend (i) incorporating more detailed vegetation types, (ii) conducting additional comparisons with specific ground data tailored to mountain environments, and (iii) performing pixel-resolution analysis to better understand the impacts of topography in future investigations.

5. Conclusions

In this study, we conducted a comparative assessment of the Plant Phenology Index (PPI) and Normalized Difference Vegetation Index (NDVI) in mountainous regions, specifically focusing on the European Alps. NDVI and PPI values were extracted from the MultiSpectral Instrument (MSI) aboard the Sentinel-2 satellite constellation for the years 2018 to 2022.

Regarding seasonal dynamics, NDVI exhibited notable signal noise during the snow seasons (late autumn to early spring) and for coniferous trees. In contrast, PPI demonstrated lower noise levels.

In terms of seasonal parameter retrieval, we observed close concordance between the two indices, with minor discrepancies in start and end-of-season estimation for deciduous trees, grasslands, and shrublands. However, substantial differences were noted for peak-of-season estimation and all seasonal parameters pertaining to coniferous trees.

Our findings align with existing literature and offer novel insights into vegetation monitoring in complex environments. Overall, this study’s primary contribution lies in exploring the potential of PPI for vegetation dynamics monitoring in mountainous regions of the European Alps, providing a robust alternative to address NDVI limitations in these environments characterized by seasonal snow and coniferous forests. As such, the PPI represents a good candidate for linking long-term regional vegetation trends to year-to-year phenological and growth responses to climate change across diverse topoclimatic contexts and life forms in mountainous environments.