Earth Observation for Phenological Metrics (EO4PM): Temporal Discriminant to Characterize Forest Ecosystems

: The study of vegetation phenology has great relevance in many ﬁelds since the importance of knowing timing and shifts in periodic plant life cycle events to face the consequences of global changes in issues such as crop production, forest management, ecosystem disturbances, and human health. The availability of high spatial resolution and dense revisit time satellite observations, such as Sentinel-2 satellites, allows high resolution phenological metrics to be estimated, able to provide key information from time series and to discriminate vegetation typologies. This paper presents an automated and transferable procedure that combines validated methodologies based on local curve ﬁtting and local derivatives to exploit full satellite Earth observation time series to produce information about plant phenology. Multivariate statistical analysis is performed for the purpose of demonstrating the capacity of the generated smoothed vegetation curve, temporal statistics, and phenological metrics to serve as temporal discriminants to detect forest ecosystems processes responses to environmental gradients. The results show smoothed vegetation curve and temporal statistics able to highlight seasonal gradient and leaf type characteristics to discriminate forest types, with additional information about forest and leaf productivity provided by temporal statistics analysis. Furthermore, temporal, altitudinal, and latitudinal gradients are obtained from phenological metrics analysis, which also allows to associate temporal gradient with speciﬁc phenophases that support forest types distinction. This study highlights the importance of integrated data and methodologies to support the processes of vegetation recognition and monitoring activities.


Introduction
Since the 1990s, under the boost of phenology networks establishment and monitoring programs initiatives, the interest in phenology has significantly increased [1][2][3][4]. Plant phenology studies the succession of plant cycle phases (e.g., leaf unfolding, flowering, blooming, leaf fall) and their development in relation to abiotic and biotic factors [5] especially with regard to meteorological components such as temperature, rainfall, humidity, radiation, exposure [6,7]. The role of phenology as a bioindicator of climate change, has been acknowledged by the World Meteorological Organization (WMO) and the Intergovernmental Panel on Climate Change (IPPC) [6,8], and the shifts in the timing of phenological stages due to climate changes have been widely recognized [9][10][11][12][13][14][15]. Phenology tracks annually recurring events in ecosystems such as plant germination, flowering, growth, and fruit maturation. Since these events are triggered predominantly by temperature, phenology has emerged as a key asset in identifying current fingerprints of climate change in nature, especially since recent warming is mirrored by significantly advancing spring events [15]. At the same time, the seasonal life cycle events throughout the year have great relevance in many fields, such as biodiversity (e.g., detection of ecosystem disturbances), cultivated species of agricultural interest (e.g., timing of management activities), forestry (e.g., monitoring forest health), and in relation to human health (e.g., providing pollen information) [2,16,17].
Nowadays, it is possible to identify time shifts and seasonal changes thank to a large amount of data set available, in some cases time series of phenological observations date back to the 19th century [18] and more [19]. This has fostered the establishment of phenological networks throughout the world, mainly based on field observations of standard growth stages of plant cycles [20], such as budburst or leaf senescence. However, ground phenology detection has some shortcomings as it is labor intensive, it involves limited geographic extents, and a limited number of species investigated [21][22][23]. On the other hand, major advances in near-surface remote sensing and satellite remote sensing have supported the developments of new methodologies of data collection and analysis from a regional to a continental scale [21,24,25]. Especially satellite remote sensing increases the range of the scale of analysis moving from plant species phenological phases to the detection of land surface phenology at the landscape scale [26,27]. It enables to investigate the spatio-temporal patterns of plant phenology and their relationship with environmental variability and climatic drivers.
The detection of phenological events by satellite remote sensing needs a dense time series of data for the purpose of being able to ascertain rapid phenological changes and to overcome the problem of cloud covered images [25,28]. To this end they made a strong contribution the Advanced Very High Resolution Radiometer (AVHHR) and Moderate Resolution Imaging Spectrometer (MODIS/Terra) sensors, providing a global coverage of very high temporal (twice daily and daily, respectively) and coarse-to-moderate spatial (1 km and 250 m, respectively) resolution observations for more than 20 years and over large geographic regions [29][30][31][32]. At present, Sentinel-2 satellites provide global acquisitions of high resolution (10-20 m) and high-revisit frequency (5 days) multispectral images enabling more detailed-scale phenological studies [33,34].
Phenology can be investigated by exploiting the optical response of vegetation through the seasonal variations of spectral and biophysical indices [35][36][37] by analyzing the vegetation index time series. Vegetation indices are the basic information for phenological researches, e.g., by analyzing the vegetation index time series to derive information [38][39][40][41][42], by assessing temporal statistics [43,44], or by estimating phenological metrics (as the onset of greenup, the length of the growing season, and the offset of the season) [45][46][47].
There are various methods to estimate the phenological metrics from vegetation index time series, namely threshold-based methods, smoothing functions, empirical equations, rate of change, change detection, phenology matching, simulation-based methods, and derivatives of vegetation greenness curves [22,25,[48][49][50][51]. Data smoothing plays a key role in the estimation of phenological metrics, since it simplifies vegetation index time series curves using empirical methods, curve fitting method or data transformation methods [25]. The successful retrieval of phenological metrics depends on the availability of robust algorithms that are capable of processing vegetation time series while minimizing atmospheric noise and sensor-related errors [52]. The aim is to extract key information from the curve of phenological time series, such as the slope of the curve or the length of the plateau, and link them to a specific phenological response.
In view of climate change mitigation and adaptation measures to be adopted by different countries, high resolution phenology trends will contribute to improve a number of monitoring activities, such as ecological status, environmental conditions, climate change impacts on ecosystems, detailed land-use/land-cover mapping, the estimation of carbon storage, vegetation responses to disturbances, cropland practices, and urban ecosystem assessments [53]. Along with the availability of high resolution satellite time series and proven methodologies to extract temporal information from satellite acquisitions, comes the need for procedures to generate high resolution Earth observation derived phenological metrics that could serve a wide range of applications. In the last decade, Sentinel-2 MSI satellite sensors have contributed significantly in encouraging vegetation investigations which have led to develop methods and products [34,54,55]. In particular, the availability of high spatial resolution and dense revisit time satellite observations, opened the way to high resolution phenological metrics estimation, representing a promising tool for the study of terrestrial ecosystems and species-specific phenology. Phenological metrics are able to synthetize key information from time series and can be used to discriminate forest typologies on the basis of the observation and measuring of the ecosystem processes responses to environmental gradients.
An automated and transferable procedure that combines available validated methodologies to exploit full satellite Earth observation time series to produce rapid and updated information about plant phenology, without any a priori information and without an excessive simplification of the temporal curve, is here presented. The aims of this study are: (i) to present a procedure to estimate phenological metrics and temporal statistics from high spatial resolution satellite images, exploiting a dense and smoothed time series of the Leaf Area Index (LAI); (ii) to investigate if smoothed vegetation curve, temporal statistics, and phenological metrics may contribute as temporal discriminants to characterize forest ecosystems in Italy.

Study Area
Italy is situated in the Mediterranean basin and covers about 300,000 km 2 ( Figure 1). It is characterized by two main mountain ranges (the Alps and the Apennines), hilly zones, river valleys, and a coastline of about 7600 km. The physiography and geographical position have led to it having heterogeneous features in terms of climate (ranging from a Mediterranean climate in the south and along the coastline, to a Temperate climate in the north and in the mountains) and in terms of land use along latitudinal and elevation gradients [56,57].
Land use is mainly represented by agricultural (52%), natural (43%) and artificial (5%) areas. Most of the natural areas are forests (28%) which are diversified into multiple habitat types.

Forest Habitat Types
The European Vegetation Archive (EVA) dataset [58,59] was used in this work. The archive includes georeferenced vegetation plots with species list and cover abundance. The sub-dataset used in this study takes into account vegetation plots smaller than 200 m 2 and with a geographic localization accuracy ranging from 0.1 to 30 m. Following the EUNIS (European Nature Information System) II and III level hierarchical classification nomenclature [60] 14,385 plots from the EVA dataset were selected and classified (Table 1) [61]. Figure 1 shows the spatial distribution of forest stands along the Italian peninsula at II level EUNIS classification. These forest types at EUNIS level II were used to estimate the smoothed vegetation curve and the temporal statistics, while regarding the phenological metrics estimate, only deciduous broadleaved forest types (T1) at EUNIS III level were considered from the vegetation plot records of the archive, since for the evergreen forest types is more difficult to identify the timing of the onset of greenness and senescence due to the lower temporal fluctuation of the vegetation index time series [61,62] and snow prevalence during wintertime in mountain ecosystems [55], hiding evergreen forests.

Satellite Data
Sentinel-2 (S2) satellite acquisitions acquired during the year 2019, with cloud cover lower than 90%, were gathered for the study area (around 7000 images, 61 granules, 9 Terabyte of data stored). The high spatial resolution (10 m, 20 m and 60 m), the high revisit time (5 days with two satellites), and the 13 spectral bands (from the visible to shortwave infrared) are the characteristics of the S2 Multi-Spectral Instrument (MSI) sensor. The images, distributed by Theia (MUSCATE format) as the bottom of the atmosphere (BOA) reflectance, orthorectified, terrain-flattened and atmospherically corrected with MACCS-ATCOR Joint Algorithm (MAJA) [63,64], were processed for spatial resampling of the spectral bands at 20 m and masked for invalid pixels (cloud, cloud_cirrus, cloud_shadow, topographic_shadow, snow, edge, and sun_too_low). Figure 2 shows the overall data processing workflow. The LAI biophysical index, defined as half of the total green leaf area per unit ground surface area, was calculated from the S2 images using the biophysical processor [65,66] available in SNAP software version 7. From the LAI time series temporal variables were calculated, spatially co-registered using the AROSICS algorithm [67], in order to deal with weak spatial coherence of S2 time series processed using processing baselines 01.xx and 02.xx [68], and stacked in a large multidimensional datacube [68]. Then the LAI time series were temporally smoothed and daily interpolated. Complementary temporally explicit information defined "weights", were generated in order to account for residual sub-pixel cloud contamination [69] and later updated during the data processing. weights values were derived from MSI B2 spectral band centered at 492 nm, and assigned following the rules: weights = 0.1 for B2 > 0.25; weights = 0.5 for B2 > 0.18 and B2 ≤ 0.25; weights = 1.0 for B2 ≤ 0.18. Figure 2. Flowchart of the main phases of the procedure to generate phenological metrics datasets and to perform multivariate statistical analysis. White ellipse refers to input data, grey boxes refer to data processing phases, black ellipses refer to output products.

Satellite Data Processing
The temporal smoothing and gap-filling aim to produce a vegetation index time series with evenly spaced timesteps, by means of a pixel-based approach consisting of four steps: (i) small drops removal; (ii) daily interpolation; (iii) weighted second-order polynomial local fitting; and (iv) Whittaker smoother. First, small drops in the temporal series were removed using a two-pass moving window filter. Daily interpolation was then performed with the R cran "stinepack" package [70], that use the Stineman algorithm [71] a method leading to much less tendency for "spurious" oscillations than traditional interpolation methods based on polynomials, such as splines. weights information was then updated with value "0.1", assigned to interpolated missing values. Later, a weighted second-order polynomial local fitting using a moving window ranging over 7 non-missing temporal observations was used to apply a fine smoothing to the time series. The determination of weights for each time step in the time series, like the application of Savitsky-Golay filter, allows to produce an upper vegetation index envelope, reducing the small drops [72].
A first-order weighted Whittaker smoother [73,74] was finally applied in order to smooth the vegetation time series; in particular, the initial and final time series observation falling outside polynomial fitting local window, using the R cran 'ptw' package [75]. Lambda value used for the smoothing was set to 10 days, after performing a sensitivity analysis to tune this value and produce consistent results when smoothing time series with a different S2 revisit time of 2-3 days (multiple orbit acquisitions on few tiles) or 5-10 days (not shown).
After the smoothing procedure, the LAI daily time series were used to estimate the phenological metrics (Table 2) using the method described in [76]: the annual temporal statistics, specifically average (LAI_avg), minimum (LAI_min), maximum (LAI_max), delta (LAI_delta), standard deviation (LAI_std), and the seasonal temporal statistics, specifically winter average (LAI_djf_avg), winter minimum (LAI_djf_min), winter maximum (LAI_djf_max), summer average (LAI_jja_avg), summer minimum (LAI_jja_min), summer maximum (LAI_jja_max). Table 2. List and description of the estimated phenological metrics used. Time field reports the units used for the estimation (x-axis in Figure 2); Value field specify if vegetation index is computed for the corresponding phenological metric (y axis in Figure 2 Plant phenology is characterized by phenophases that follow each other, from dormancy (or sowing in agriculture practices) through the start of season (characterized by the onset of greening with leaf unfolding), peak (or heading), maturity phase, flowering, senescence (characterized by curtailing of chlorophyll production that reveals various accessory pigments), leaf fall, and finally, dormancy.
The method to estimate phenological metrics is based on a combination of local maxima in the first derivative [22,76] (Figure 3). It first identifies the vegetation peak, corresponding to absolute maximum time series value. Peak is used to initialize the identification of maximum increase (temporarily happening before the peak) and maximum decrease (temporarily happening after the peak) of time series first derivative. Maximum greenup rate (hereafter greenup) and maximum senescence rate (hereafter senescence) are defined as slopes of recovery and senescence lines tangent to the time series. The intersections among these lines and baseline and maxline identify the four phenophases in the original formulation (UD: upturn date; SD: stabilization date; DD: downturn date; RD: recession date) [76]. Starting from these phenophases, phenological metrics were extracted ( Figure 3). Dates were expressed as absolute number of days (i.e., DoS and LMP phenological metrics) or as both calendar date and day of year (DoY).

Phenological Metrics Accuracy Assessment
In order to verify the accuracy of the estimated phenological metrics with ground phenology observations [51,77], the resulting S2 estimates were compared with the digital images of the PhenoCam network observations (https://phenocam.sr.unh.edu/webcam/, accessed on 11 March 2021) available for natural ecosystems in the study area ( Figure 1, Table 3). The PhenoCam Dataset V2.0 is a fully processed, quality-controlled, and curated data set which is made freely available through the ORNL DAAC [78,79]. The daily 90th percentile of Green Chromatic Coordinate (GCC 90 ), a vegetation index derived from Phe-noCam photographic images that quantifies the greenness relative to the total brightness, has been used as reference time series and processed with the same procedure used to analyze S2 time series, starting from temporal smoothing processing step. Accuracy of the estimated phenological metrics was assesses by computing the Mean Error (ME), the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE) and the correlation coefficient (r) between reference field data and S2 estimates. Since phenological metrics were derived from different time series (i.e., GCC 90 and LAI), with a nonlinear relationship, only date estimates (expressed as DoY or absolute number of days) were compared for all the available PhenoCam time series in the period 2016-2020. S2 acquisitions available in the period 2016-2020 were collected for location correspondent to PhenoCam sites, in order to generate multi-year time series for the comparison with ground observations.

Multivariate Analysis
Discriminant Function Analysis (DFA), Canonical Correlation Analysis (CCA) and Linear Discriminant Analysis (LDA) were selected to analyze the contribution of the three generated datasets (the smoothed vegetation curve, the temporal statistics, and the phenological metrics) to the characterization of forest habitat types ( Table 1). The DFA is a statistical method that separates objects into classes by performing a multivariate analysis to highlight the differences among groups, namely the categorical response variable, and the variables needed to describe these differences, namely the predictors [80]. Four different forest habitat classes were considered for DFA analysis, corresponding to the broadleaved deciduous (T1), broadleaved evergreen (T2), needleleaved evergreen (T3) and needleleaved deciduous (T34) plant functional types. The CCA was used to explore the relationship between two sets of variables combining them in order to extract the variance of the data matrix by a limited number of independent axes (i.e., "canonical roots", or "roots") [81]. Linear Discriminant Analysis (LDA) is a dimensionality reduction technique, to reduce the number of dimensions in a dataset while retaining as much information as possible, providing the best possible separation between the groups [82].
The georeferenced points of vegetation plots were used to perform a spatial query over the variables generated from satellite time series, and precisely all the plots (14,385) over the smoothed vegetation curve and temporal statistics, and the T1 plots (8328) over the phenological metrics.
The data analysis was performed according to the following steps: 1.
Regarding the smoothed vegetation curve and the temporal statistics variables (predictors), a DFA was executed to identify the variables able to discriminate the forest types (response variables) at the EUNIS II level; 2.
As for the phenological metrics, a two steps analysis was performed for T1 EUNIS III level classes. Since the presence of two sets of variables, firstly a CCA was run on the phenological metrics dataset illustrated in Table 2 and constituted of 10 variables of LAI values and 9 variables of date values. The resulting independent and representative axes were then used to perform a LDA to discriminate the deciduous broadleaved forest types.
The analyses were performed using GDAL libraries, R programming language and libraries, QuantumGIS, SeNtinel Application Platform (SNAP), NetCDF Operators, and Climate Data Operators software. Results of the statistical analysis, conducted to compare the phenological metrics estimated from S2 time series with reference field data estimated from PhenoCam time series, located in the Alpine region, are reported in Table 4. Figure 5 shows the distribution of phenological metrics dates estimated from ground observation and from S2 time series.

Results
Two EGS dates clusters can be observed in Figure 5, likely related to the differences between tree species and grasslands, the latter characterized by an early yellowing phase. PoS DoY estimated from S2 generally occurs later than the PhenoCam estimated one. This may be related to higher GCC 90 signal saturation than LAI, and result in higher error metrics (Table 4) as compared to other assessed phenological metrics. Error metrics for DoS and LMP, estimated from a pair of phenological metrics estimates, may be higher than others due to the uncertainties in both greenup and senescence phenophases, that can increase error. Figure 6 shows the results of multivariate analysis. Regarding the smoothed vegetation curve, the first two factors (F1 and F2) explain the 91.4% of the total variance (Figure 6a) for the EUNIS II level forest types and the 44.3% of the total variance for the EUNIS III level forest types (Figure 6b). In Figure 6a, the F1 (72.6% of the total variance) reveals a seasonal gradient with the deciduous forest types placed in the negative side of F1 and the evergreen in the positive one. The F2 instead (18.9% of the total variance) shows a differentiation of leaf types with the broadleaved located in the negative side of F2 and the needleleaved in the positive one. Figure 6b shows a similar separation between deciduous and evergreen alongside the F1, and broadleaved and needleleaved alongside F2, except for T33 (Mediterranean mountain Abies forest) that is a mixed forest type.   As for the temporal statistics, F1 and F2 explain the 99.6% of the total variance (Figure 6c) for the EUNIS II level forest types and the 87.5% of the total variance for the EUNIS III level forest types (Figure 6d). The distribution of the forest types with respect to the four quadrants of the DFA biplot is in line with what was reported for the LAI time series. Indeed, it is possible to recognize the same gradients in Figure 6c,d. Differently from the LAI time series biplots, that does not clearly express evident gradients, it is possible to find the contribution of the temporal statistics to the forest types discrimination. In Figure 6c, the F1 (94.4% of the total variance) reveals the summer productivity (LAI_max, LAI_jja, and LAI_std) and discriminates the deciduous forest types with T1 and T34 located in the negative side of the F1. T2, which represents the forest type with the maximum availability in terms of number of days of photosynthetic surface, is instead discriminated by LAI_avg and LAI_djf_max. The biplot of Figure 6d indicates a leaf-type productivity with T2 forest types located in the positive side of F1 (75.6% of the total variance) and in the negative side there are T3 and T1 forest types.
Regarding the phenological metrics, the output of the CCA consists of two sets of independent axes: SET1 axes for time values (DoY) and SET2 axes for LAI values. The results of the LDA are shown in Figure 6e,f. LD1 and LD2 explain the 81.5% of the total variance (Figure 6e), whereas LD2 and LD3 explain the 57.7% of the total variance (Figure 6f). In Figure 6e, LD1 (48.5% of the total variance) shows a temporal gradient with SET1_SCORE1, SET1_SCORE2, and SET1_SCORE3 that discriminate forest types characterized by a shorter growing season located in the right side of LD1. LD2 (33% of the total variance) is an expression of LAI values and discriminates the altitudinal gradient of broadleaved deciduous forests with beech and alder forest at the bottom of the biplot, and beech and oak at the top. Moreover, mesophilous forest types are located on the top side of the biplot. In Figure 6f, LD1 (48.5% of the total variance) shows a similar temporal gradient as in Figure 6e, whereas LD2 (9.2% of the total variance) indicates a latitudinal gradient with the alpine forest types located on the top side of the biplot (T18, T1C, and T15). To sum up, a smoothed vegetation curve accounts for seasonal gradient and leaf type discrimination, whereas temporal statistics results, although showing the same gradients, enhance forest-type discrimination, highlighting forest productivity and leaf productivity features. Phenological metrics outcomes allow to enrich the information provided by the previous two analyses by associating a temporal gradient with specific phenophases that support forest-type distinction. Furthermore, altitudinal and latitudinal gradients also arise.

Discussion
Development of Earth observation data analytics allows to analyze and systematically extract diverse sets of information from a variety of large datasets, including those observing and measuring the response of environmental and ecosystem processes. Earth observation derived phenological metrics represent a promising tool to monitor ecological status, environmental conditions, climate change impacts on ecosystems, cropland practices, and more accurately forecast crop yields. Trends of phenological shifts in both spatial and temporal scales, with consequent impact on ecosystem functioning, could be identified using high resolution satellite derived phenological metrics.
Land-surface phenology has been estimated in the last decades from medium spatial resolution high revisit time satellite observations, that allow observing on regional to global scales but have a limited representativeness for phenological changes at ecosystem or species-level. S2 satellites, with a high spatial resolution and 5 days revisit time, opened the way to high resolution phenological metrics estimation, and represents a promising tool for a variety of ecological analyses, including for example the study of terrestrial ecosystems and species-specific phenology [22]. The LAI vegetation index estimated from satellite observations has been selected as source information to derive phenological metrics since it represents a biophysical parameter, and it is less affected by signal saturability in areas with high vegetation coverage. Alternatively, NIRv [83] and kNDVI [84] vegetation indices can be considered since they have low signal saturability in comparison to commonly used NDVI and EVI vegetation indices.
Retrieving plant phenology from time series of satellite Earth observation vegetation indices has been widely investigated and applied in the past two decades. USGS EROS and Copernicus initiatives [54,85] contributed to the development of operational products for landsurface phenology monitoring. TIMESAT, an algorithm implementing least-squares methods for processing time series of Earth Observation data, has been largely used to estimate phenophases. The implemented Savitzky-Golay filtering works very well with time series relatively unaffected by noise, while the asymmetric Gaussians, classified as semi-local method, force the data into the fixed functional form and they are able to follow also more complex behaviors, such as a rapid increase followed by a decreasing plateau [86]. The asymmetric Gaussian method has been found to be less sensitive to noise and to give better predictions for the beginnings and ends of the seasons [86]. The enhanced TIMESAT algorithm, using the third derivative, is relative stable in determining greenup and senescence, no matter whether vegetation index is changing quicker or slower [32].
Among the existing methodologies, the ones used in this study to temporally smooth time series and to estimate phenology metrics have been selected considering requirements and previous works in literature. The requirement for a generalized procedure, that can be applied to other vegetation indices and other satellite sensors, is the use of methods without any a priori information. Keeping smoothed temporal signal as close as possible to actual observations requires the use of a local curve fitting, to avoid an excessive simplification of the temporal curve methodologies and to allow less altered estimates for some of the phenological metrics (e.g., plateau_slope and STI). Finally, the use of a co-registration method to improve high spatial resolution satellite observations is advisable. The proposed procedure uses local curve fitting and local derivatives to identify phenophases, operating without thresholds or a priori information. Signal filtering is a very important processing step because it makes the vegetation time series interpretable to retrieve phenological infor-mation. The methodology temporally interpolates S2 time series in order to homogenize seasonal trajectories of vegetation indices. All the available satellite observations need to be optimally used, even if heavily hampered by clouds. In this study, the removal of invalid pixels (e.g., clouds, topographic shadows) prior the analysis assumes a key role in reducing factors affecting the quality of observations, that typically results in reducing vegetation indices values. The selection of S2 dataset with rigorous cloud detection algorithm [87], the use of weighted smoothing procedure and the temporal image co-registration enhance vegetation indices time series integrity.
Phenophases are estimated using the method presented in [76] and implemented in Phenopix R package, which offers a suite of standardized and reproducible processing code, designed for the extraction of phenological information from time-lapse digital photography of vegetation cover [22]. The method allows the extraction of greenup and senescence maximum rate date from time series smoothed using local curve fitting, diversely from other methodologies based on derivatives of the vegetation index seasonal trajectory that require extremely smoothed time series. The proposed approach distinguishes between Start of Season (SoS)-End of Season (EoS) and Start of Growing Season (SGS)-End of Growing Season (EGS) metrics, as compared to similar approaches. The start and the end of growing season is intended to represent seasonal bounds of photosynthetic tissues development phases. Diversely, the end of season comes at the end of the period of leaf-mineral nutrient remobilization during leaf senescence, when the plant is curtailing chlorophyll production revealing various accessory pigments, determining yellowing or browning of leaves, and consequently, driving the starting of fall foliage [61]. Vegetation normally changes more quickly during greenup than senescence [32], as noticeable from Figure 4. SoS phenological metric is of particular interest in the agriculture field since it represents the timing when vegetation index is at minimum before the onset of greenness. It should be related to pre-sowing plowing, a soil management practice that deals with soil preparation for a new crop, also utilized as post-harvest weed infestation control.
In forest ecosystems, minimum LAI is likely dependent on both evergreen overstory and/or understory species within a mixed-forest pixel [35]. The effect of deciduous understory species may affect surface reflectance very early in the growing season before growth of overstory canopy occurs. As a result, an early detection of onset of greenness for tree species within a mixed-forest pixel can be estimated [35]. To reduce such effect, amplitude metric is computed considering the difference between vegetation peak and vegetation baseline, in order to eliminate the effect of evergreen species on the minimum vegetation index value.
The association between the estimated phenological metrics dates and ground observations was evaluated. To assess the accuracy of the phenological retrievals, ground-truth information on phenological transition dates would be required [88]. Although subjective ocular estimates were shown to have a close match to remote sensing indices, they are time-consuming as they require the frequent presence of an observer in the field [89]. To overcome this, the use of digital cameras that automatically take several pictures per day, with a strong focus on forest ecosystems [21], allows to collect time series of greenness chromatic coordinates, that can then be subjected to similar phenology extraction methods as for satellite imagery to identify phenological transition dates [22]. Validating phenology metrics derived from satellite data product may be difficult due to vegetation heterogeneity [32]. Data quality and phenology retrieval methods can be source of uncertainties in phenological metrics estimates [90], and its accuracy remains to be validated globally. Error metrics calculated using ground observations resulted in a MAE of about 15 days for the various metrics analyzed. This is consistent with accuracy reported in other studies [34,91] and can be caused by many error sources. For example, the use of different indices to estimate phenological metrics, LAI and GCC 90 with the latter more affected by signal saturability (Figure A1), may have contributed to generate differences in the estimated values. In fact, previous research studies demonstrated that phenological metrics derived from NDVI and EVI for the same geographic area are not in 100% agreement with each other [32].
Forest plant communities, together with their own typical floristic composition, show exclusive phenological dynamics recognizable by vegetation indices time series [42]. Both vegetation indices time series [42,61] and land surface phenological metrics [61,92] has been successfully used to classify plant communities and natural habitats, discriminating vegetation patterns dominated by species with similar phenology features. Multivariate statistical analysis conducted in this research study allowed to understand how the three generated datasets (the smoothed vegetation curve, the temporal statistics, and the phenological metrics) have the capacity to serve as temporal discriminants to detect forest ecosystems processes responses to environmental gradients. Results from multivariate analysis demonstrate how temporal statistics and phenological metrics are representative of the time-related variability, can synthetize key information from satellite time series, reducing data dimensionality, and thus can be used as temporal discriminants for forest ecosystems classification and mapping. Specifically, smoothed vegetation curve and temporal statistics are able to highlight seasonal gradient and leaf-type characteristics to discriminate forest types, with additional information about forest and leaf productivity provided by temporal statistics analysis. Furthermore, temporal, altitudinal, and latitudinal gradients are obtained from phenological metrics analysis, which also allows to associate temporal gradient with specific phenophases. High spatial resolution smoothed time series and phenological metrics open up to the provision of novel temporal information about forest phenology anomalies and useful monitoring system to scrutinize spatio-temporal patterns of forest disturbance [62].
In the frame of ecosystem monitoring (e.g., conservation status), and specifically for grassland management in agricultural areas (e.g., fodder), the systematic capturing of cutting times would be highly relevant for the surveillance of areas granted by EU funding programmes [93], under policies on rural development through funding and actions (EU's Common Agricultural Policy-PAC). The capacity of mapping crop types using high spatial resolution phenological metrics estimated using the proposed procedure has been demonstrated in [94] for heterogeneous, small, and fragmented agricultural systems. In the context of crop-type mapping and the monitoring of agricultural practices, synthesizing information to fewer phenological metrics would facilitate image data processing, for example, image segmentation by reducing time series dimensionality.
While inter-annual variability of phenological metrics can be evaluated even at a local scale, and analysis on continental scales can detect spatial variability in phenology across climate gradients. In fact, vegetation phenology is highly sensitive to climate conditions and is a climate change fingerprint [15]. Warm and cold spells, which are not single extreme events but can be regarded as a compound extreme, such as a persistence of weather conditions, have impacts on the onset of phenological seasons, that differed strongly depending on species, phase and timing of the event [15]. Phenology also controls many feedbacks of vegetation to the climate system by influencing the seasonality of albedo, surface roughness length, canopy conductance, and fluxes of water, energy, CO2 and biogenic volatile organic compounds [95]. Knowledge of how geo-morphological and bio-climatical conditions affect phenological behavior is valuable information to model the potential effects of climate change and to estimate the future adaptation of plant growing in different geographical areas.
The proposed approach has been demonstrated using S2 LAI time series. Other vegetation indices and other satellite sensors or virtual constellation of sensors can be used to estimate phenological metrics with the presented procedure. New perspectives concerning monitoring of plant phenology can benefit from dense time series generated from virtual constellations, such as the Harmonized Landsat and Sentinel-2 (HLS) dataset initiative, or synergies with Synthetic Aperture Radar (SAR) sensors satellite observations. SAR time series have a strong seasonal signal in VH radar backscattering coefficient [96], as well as the ratio between VV and VH coefficients, that allows the monitoring of phenology in deciduous forests independently of the cloud cover [97].

Conclusions
Earth observation derived phenological metrics to represent a promising tool to monitor ecological status, environmental conditions, climate change impacts on ecosystems, cropland practices, and more accurately forecast crop yields. The proposed automated and transferable procedure combines available validated methodologies to exploit full satellite Earth observation time series without any a priori information, to produce rapid and updated information about plant phenology. Estimated phenological metrics have been validated using in situ PhenoCam observations with satisfactory results.
This study shows that integrated data and methodologies based on plant phenology may be an effective tool to generate Earth Observation derived temporal discriminants, very useful to characterize forest ecosystems, and may help the processes of vegetation recognition and classification, in addition to support monitoring activities of natural ecosystems and agro-forestry systems.