Estimation of Gross Primary Productivity (GPP) Phenology of a Short-Rotation Plantation Using Remotely Sensed Indices Derived from Sentinel-2 Images

: This study aimed to understand which vegetation indices (VIs) are an ideal proxy for describing phenology and interannual variability of Gross Primary Productivity (GPP) in short-rotation coppice (SRC) plantations. Canopy structure- and chlorophyll-sensitive VIs derived from Sentinel-2 images were used to estimate the start and end of the growing season (SOS and EOS, respectively) during the period 2016–2018, for an SRC poplar ( Populus spp.) plantation in Lochristi (Belgium). Three different filtering methods (Savitzky–Golay (SavGol), polynomial (Polyfit) and Harmonic Analysis of Time Series (HANTS)) and five SOS- and EOS threshold methods (first derivative function, 10% and 20% percentages and 10% and 20% percentiles) were applied to identify the optimal methods for the determination of phenophases. Our results showed that the MEdium Resolution Imaging Spectrometer (MERIS) Terrestrial Chlorophyll Index (MTCI) had the best fit with GPP phenology, as derived from eddy covariance measurements, in identifying SOS-and EOS-dates. For SOS, the performance was only slightly better than for several other indices, whereas for EOS, MTCI performed markedly better. The relationship between SOS/EOS derived from GPP and VIs varied interannually. MTCI described best the seasonal pattern of the SRC plantation’s GPP (R 2 = 0.52 when combining all three years). However, during the extreme dry year 2018, the Chlorophyll Red Edge Index performed slightly better in reproducing growing season GPP variability than MTCI (R 2 = 0.59; R 2 = 0.49, respectively). Regarding smoothing functions, Polyfit and HANTS methods showed the best (and very similar) performances. We further found that defining SOS as the date at which the 10% or 20% percentile occurred, yielded the best agreement between the VIs and the GPP; while for EOS the dates of the 10% percentile threshold came out as the best.


Introduction
Short-rotation coppice (SRC) plantations provide multiple economic and ecological benefits [1]. They are ideal for producing biomass for energy [2] and can, therefore, substitute for fossil fuels to reduce anthropogenic CO2 emissions [3]. They enrich farm-scale biodiversity [4], increase soil organic carbon stock [5,6], improve groundwater quality [7] and help in preventing pests and diseases. Over 16% of woodlands in Europe were managed as coppice in 2000, covering a global area of about 23 million ha [8]. In the last 20 years, SRC plantations increased globally to cover an area of 31.8 million ha in 2016 [9]. Poplars (Populus spp.) are widely used as bioenergy plantations and for wood production in the Northern Hemisphere representing 99% (31.4 million ha) of SRC area worldwide.
Even though the important role of SRC plantations in bioenergy productions and carbon sequestration is well recognized by the scientific community [10,11], data on SRC productivity (e.g., Gross Primary Productivity (GPP)) and phenology are scarce and limited to small experimental areas.
The ecosystem phenology, also defined as carbon flux phenology, describes the seasonality of the gross photosynthesis (i.e., GPP) of an ecosystem. The starting day of the photosynthetic season (SOS) occurs when an ecosystem turns from a carbon source into a carbon sink in spring. Contrary, the end of the season (EOS) identifies the decrease of photosynthetic activity at the beginning of senescence. A shift in the SOS and EOS of a growing season regulates annual GPP [12,13]. A longer growing season increases the annual productivity of deciduous forests [14,15]. The timing of bud bursts and sets depend on local environmental conditions. The plant growth strategy of deciduous trees in temperate regions exploits maximally optimal spring growth conditions [16]. SRC plantations allocate carbon to different tree components at the beginning of the growing season [17] as reported also for beech in Campioli et al. (2011). The photosynthetic activity shows a substantial decline by early mid-September in temperate deciduous species [18]. Contrary, high carbon uptake rates were observed until the end of September for SRC plantations [17,19] and poplar [20]. Furthermore, for agroecosystems such as crops and SRC plantations, management practices also play an important role in determining the annual production. However, few studies focused on the understanding of the effect of coppicing on SRC phenology [21,22]. A clear understanding of how phenology affects the productivity of SRC plantation can effectively support farmers in optimally managing their fields to obtain the most profit from crop production. In addition, it is necessary to improve the model parameterization of global productivity models.
In situ observations provide a robust measurement of plant development and phenology [22], but they can be laborious and expensive when a lot of sites need to be sampled. The eddy covariance (EC) method [23] is suitable for monitoring ecosystem productivity, aspects of phenology and their changes related to abiotic stresses (e.g., drought, management) [24,25]. However, GPP data derived from the EC method are representative of only a small area, roughly speaking 10 to 30 times the measurement height on average, around and mostly upwind of the flux tower [26] and the technique has rarely been applied at SRC plantations. Both EC and in situ observations thus have limitations when data sampling is required at large spatial scales.
Large-scale and continuous observations taken by satellite sensors make remote sensing approaches attractive alternatives for phenology assessment [27,28]. Spectral vegetation indices (VIs) can be used to derive many biophysical and biochemical vegetation parameters (e.g., chlorophyll content, Leaf Area Index (LAI)) and GPP [29][30][31][32][33][34][35]. The Normalized Difference Vegetation Index (NDVI) [36] and the Enhanced Vegetation Index (EVI) [37] (referring to Table 1) are the most commonly used VIs for monitoring canopy greenness, development and phenology [38]. Due to its long time-series, NDVI has been widely used as a proxy for GPP [39,40] and the analysis of ecosystem phenology [28,[41][42][43][44]. Research, however, has shown that both NDVI and EVI can fail to predict canopy phenology and productivity when LAI is high and in specific periods of the year, especially when canopy greenness and physiology are decoupled, e.g., during drought [45,46] and senescence [47,48]. Chlorophyll-sensitive VIs appear to be better suited for monitoring the phenology of plants and crops [49], because they are indirect proxies of canopy biochemistry. For example, green NDVI (gNDVI) considers the green band, which is more closely related to chlorophyll (see Table 1) than the red band [50], which is used in NDVI. gNDVI is better correlated with LAI than are NDVI and EVI and improves the modeling of LAI [51]. The Chlorophyll Index (ChlRedEdge) and the MERIS Terrestrial Chlorophyll Index (MTCI) also correlate with chlorophyll content [52][53][54] and should, in theory, be good proxies for phenology and GPP seasonality of crops. Pigment Specific Simple Ratio (PSSR) was reported to have a strong relationship with chlorophyll content [29,55]. Modified Chlorophyll Absorption Ratio (MCARI) is an alternative index that provides a measure of chlorophyll concentration [56,57] but is also sensitive to nonphotosynthetic components of the canopy, especially at low chlorophyll concentrations [58]. The Structure Insensitive Pigment Index (SIPI) is associated with the ratio of carotenoids to chlorophylls and allows mitigating the effects of foliar structural properties on reflectance, because it considers the NIR band [59][60][61]. Thus, a wide suite of RS-based indices can provide information on SRC phenology and possibly productivity. However, to date, the phenology and GPP seasonality of SRC plantations have not yet been widely studied via VIs, probably due to the poor spatial (1 km or larger) and spectral resolution (red and NIR mainly) of satellite sensors (e.g., MODIS/Terra or MODIS/Aqua; SPOT/VGT, NOAA/AVHRR; Landsat) [22].
A wide range of time-series smoothing methods and extracting phenological metrics (SOS, EOS) to reconstruct the seasonal patterns of noisy VI time-series have been proposed [27,62,63]. However, there is no single method that always performs better than another for smoothing VI time-series [64]. Large discrepancies in the detection of SOS and EOS have been found by comparing different retrieval methods [65][66][67] and no conclusive evidence favors one method over another. Therefore, combining multiple methods should provide better estimates of phenological metrics. Green biomass, canopy greenness and phenology [70] The recent advent of high-resolution Sentinel-2 satellite imagery, with a spatial resolution of 10, 20 and 60 m offers a new perspective in monitoring ecosystem phenology. This satellite has two main advantages over other high-resolution satellites such as Landsat. The five-day temporal resolution is an advantage for monitoring plants, and the 10 m spatial resolution is more suitable for ecological, agricultural and phenological studies [71,72]. Another interesting feature of Sentinel-2 is the availability of additional red-edge bands, which are available from other few space-borne sensors [73]. These red-edge regions play an important role in estimating the biochemical and biophysical traits of plants [74].
In this study we explored for the first time the use of eight canopy structure-and chlorophyllsensitive VIs, derived from Sentinel-2, for describing the phenology and seasonal variability of GPP, derived from EC, of an SRC poplar plantation in Belgium. The objectives of this study were threefold: (i) to understand which VI could be an ideal proxy of GPP phenology of SRC plantations, (ii) to identify optimal filtering and threshold methods for describing phenological metrics and (iii) to analyze the performance of each VI in capturing the seasonal variability of GPP.
A complete description of the experimental site and plant material has been previously published [75,76]. In brief, the plantation consists of 12 poplar (Populus) genotypes planted in monoclonal blocks at a density of 8000 trees ha-1. Poplars are planted in a double-row planting scheme with an alternating interrow spacing of 0.75 m and 1.50 m and an intrarow distance between each tree of 1.

Measurement of CO2 Fluxes and Estimation of GPP
The study site was equipped with an EC system to monitor continuously the CO2 net flux between the surface and atmosphere, also called Net Ecosystem Exchange (NEE). The system was composed of a Gill-HS50 sonic anemometer (Gill Instruments Ltd., Lymington, UK) for measuring 3D wind speed components and sonic temperature and of an LI-7200 enclosed-path infrared gas analyzer (LI-COR Inc., Lincoln, NB, USA) for measuring CO2 and H2O concentrations. Both instruments collected measurements at a frequency of 10 Hz. CO2 and H2O concentrations and fluxes of latent and sensible heat were estimated every 30 min using a standard scheme for processing the raw EC data. Briefly, flux data collected during conditions of well-developed turbulence and a good level of stationarity with sustained mechanical turbulence were used for the analysis according to CARBOEUROPE methodology [77]. The location of the tower was on the north-east edge of the plantation (Figure 1). To maximize the representativeness of flux measurements collected at the EC tower, data with wind direction between 50° and 250° were retained for the analysis [78] (Figure 1, wind rose). A more detailed description of the experimental setup, the calculations applied to the raw EC data and the criteria used to filter the entire data set have been previously reported by [11,19]. Complete details of the EC methodology, processing schemes and corrections are also available [77,79].
Other meteorological variables were continuously measured and averaged every 30 min at the site: global solar radiation (Rg) was measured using a CG3 pyranometer (Kipp and Zonen B.V., Delft, The Netherlands), air temperature and relative humidity were measured using an HMP155A thermohygrometer (Vaisala Oyj, Vantaa, Finland) and soil water content at a depth of 20 cm was measured using CS650 time-domain reflectometric probes (Campbell Scientific Inc., Logan, UT, USA).
GPP was estimated from the measured NEE by partitioning it into its two main components: GPP and ecosystem respiration ( ), following the well-established relationship = − , where Reco is modeled from NEE night-time data as a function of measured air temperature [80,81]. The full time-series of the three variables have been obtained by filling unavoidable gaps in the measured NEE time-series, using the methodology for gap-filling and partitioning described by Reichstein et al. (2005). Numerical results were obtained using the REddyProc online tool provided by the Max Planck Institute (https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWeb). We averaged half-hourly GPP values between 10:00 and 14:00 solar local time to allow for direct comparison with the Sentinel-2 satellite data acquisition. These values are indicated hereafter as GPP, while daily sums of GPP are indicated as daily GPP.

VIs
In this study, we used the European Space Agency Sentinel-2 TOC (Top Of Canopy) satellite images provided by the Terrascope platform (https://www.terrascope.be). This platform provides access to a data set of processed Sentinel-2 images covering Belgium. Sentinel-2 carries a Multispectral Instrument (MSI) that samples 13 spectral bands in the visible and NIR regions at 10, 20 or 60 m resolutions depending on the spectral band (https://sentinel.esa.int/web/sentinel/home). The whole time-series of Sentinel-2 data from 1 January 2016 to 31 December 2018 were used in this study. The long cloudy and rainy periods in this area, mainly in 2016, affected the quality of satellite time-series. After the quality check, 15 images were retained for 2016, 24 for 2017 and 23 in 2018. From the Sentinel-2 time-series an area of 8 ha corresponding to the area of the plantation was extracted first from the Sentinel-2 time-series for each available day ( Figure 1, red polygon). A total of about 800 pixels were extracted for the bands at a resolution of 10 m and 200 pixels were extracted for the bands at a resolution of 20 m (see Table 1 for more details). Secondly, using the cloud mask available for this data set, each pixel was quality checked and only cloud-free days were selected for further analysis [82]. Third, the values of the retained good pixels in the extracted area were averaged for each band, and VIs were derived. Eight VIs were calculated over all available images of the timeseries: the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI), the MERIS Terrestrial Chlorophyll Index (MTCI), the Chlorophyll Red Edge Index (ChlRedEdge), the Modified Chlorophyll Absorption Reflectance Index (MCARI), the Pigment Specific Simple Ratio (PSSR), the Structure Insensitive Pigment Index (SIPI) and the green NDVI (gNDVI). The vegetation indices were calculated with Sentinel-2 bands B02, B03, B04, B05, B06 and B08, using the formulas reported in Table 1. The Sentinel-2 bands have the following band centers; B02: 490 nm; B03: reflectance at 560 nm; B04: reflectance at 665 nm; B05: reflectance at 705 nm; B06: reflectance at 740 nm; B08: reflectance at 842 nm. B02, B04 and B08 have a 10 m spatial resolution, while B05 and B06 have a 20 m resolution.
The NDVI and the EVI are the most known and used spectral indices for monitoring canopy greenness, development and phenology [38]. However, both the NDVI and the EVI have shown to fail in predicting the canopy phenology and productivity in specific periods of the year, especially when canopy greenness and physiology are decoupled, i.e., during drought [45] and senescence [47]. For that reason, in this study, we also considered VIs that correlate with canopy pigments and, therefore, are indirect proxies of canopy biochemistry. In particular, the gNDVI is more sensitive to chlorophyll. The ChlRedEdge and the MTCI indirectly correlate with chlorophyll content [52][53][54]. The PSSR was reported to have a strong relationship with chlorophyll content [29,55] while MCARI correlates to chlorophyll concentration [56,57]; the Structure Insensitive Pigment Index (SIPI) is sensitive to the ratio of carotenoids to chlorophyll [61].

Extracting Phenology Indicators from GPP and VIs
Dates for SOS and EOS were derived from time-series for both GPP and VIs. Three smoothing methods were applied to the raw time-series data to reduce noise and fill gaps in the occasionally poor coverage: Savitzky-Golay (SavGol), Harmonic Analysis of Time Series (HANTS) and polynomial (Polyfit).
All data of the fitted time-series before the period with the highest GPP (April−June) were selected for identifying SOS while the remaining data were selected for identifying EOS. Five different threshold methods were applied to both data sets to identify SOS and EOS, respectively.

Filtering Methods
The study area was characterized by long rainy and cloudy periods, mainly in the year 2016. For that reason, we selected three smoothing functions, whose in regions with a long period of cloudiness and a high atmospheric variability are well documented [28,83,84] The SavGol filter is a robust method to reduce efficiently the cloud and snow contamination in the time-series [84,85]. It is a low-pass filter that smooths the time-series applying a polynomial regression function to data points in a moving window [86]: For all points within the time window, a polynomial function of degree was solved. In this study, the window size and the degree of the polynomial fitting function have been defined for each time-series data (i.e., VIs and GPP). In every window, a new polynomial was fitted. Subsequently, the window moved to the next part of the input dataset, a new local polynomial function was fitted to the original data in this window, and a new data point was calculated. This process was repeated for all remaining time windows. The parameters used to smooth both VIs and GPP time-series are reported in Table 2. We obtained these parameters using Python 3.6 package scipy.signal and the function SavGol filter (https://scipy.org/). HANTS filter employs the Harmonic Analysis of Time Series, adapted from the fast Fourier transform [87]. A least-square fit was applied iteratively on all data points of each time-series, the observations were compared to the fitting curve and outliers removed from the time-series [88]. This analysis was done using the HANTS script in MATLAB 6.3 (MathWorks Inc., Sherborn, MA, USA) developed by Mohammad Abouali (https://mabouali.wordpress.com/projects/harmonic-analysis-oftime-series-hants/). For both GPP and VI time-series we used the following default parameters: valid range minimum = 0.0; fit error tolerance = 5.0; degree of over determinedness = 1; delta to suppress high amplitudes = 0.1. The length of the base period was settled to 2. For GPP, the total number of the time-series and the number of frequencies to be considered above the zero frequency were fixed to the total number of days in a year (365 or 366). For VIs, these two last parameters corresponded to the number of Sentinel-2 images available for each year.
Polyfit method was used to apply a -degree polynomial function (Equation (1)) to fit the whole time-series of data without applying a moving window. The coefficients of the polynomial function were obtained by the least-square fitting method (see Table 2). This analysis was performed using Python 3.6 package scipy.signal and the function polyfit (https://scipy.org/).
The accuracy of each filtering method in describing GPP seasonality for each year was evaluated by the root mean square error (RMSE) between fitted GPP and GPP simulated by a linear model obtained between fitted GPP and VI time-series. For each filtering method and year, this analysis was applied twice for the two periods used to extract the phenological metrics by thresholds: April−June and July-December. Table 2. Parameters used for Savitzky-Golay (SavGol) and polynomial (Polyfit) methods for filtering time-series of Gross Primary Productivity (GPP) and VIs (see VIs definitions in Table 1) per each year.
-degree of polynomial function; -data points in the moving window used for the SavGol method. gNDVI 5 3 5

Threshold Definitions
The extraction of phenological metrics (SOS and EOS) from VI time-series is highly uncertain. Large discrepancies in the detection of SOS and EOS have been found by comparing different retrieval methods [65][66][67]89]. We tested three state-of-the-art thresholds (first derivative, percentile and percentage) to perceive their differences and understand which one is the ideal threshold for GPP and VIs.
First derivative. The SOS occurs when the first derivative value of the time-series curve changes from negative to positive, while for EOS is when it changes from positive to negative [90][91][92]. The max and minimum of the first derivative correspond to the faster greening and declining (senescence) rates respectively [93].
Percentile. The percentile method does not have a universal definition. The two data sets for SOS and EOS were ordered from lowest to highest and the percentile was calculated as [94]: where is the reference percentage and is the number of data points [95,96].
Percentage. The percentage considers the difference between the highest and lowest values. The data of each data set were first ordered from lowest to highest, similar to the percentile method [97]. The percentage was calculated as: where MaximumN and MinimumN are the highest and lowest values among the time-series data (e.g., GPP or VIs) and is the selected percentage (e.g., 10%, 20%). After applying the equation to both data sets, we selected the corresponding value among the data set (e.g., parameters such as GPP or VIs) equal to or larger than the result of this equation. The days of the year associated with this value were selected as SOS and EOS for that time-series (e.g., GPP or VIs).
For both percentile and percentage threshold methods, the reference percentage (in Equations (2) and (3)) indicated the point in time when a defined fraction of the seasonal amplitude has been reached. We set to 10% and 20% because thresholds between 10% to 25% provide more accurate estimates of phenological metrics as previous works have demonstrated [98,99].

Statistical Analysis
To define which smoothing function or threshold of VIs was most suitable for detecting SOS and EOS of GPP, two statistical analyses were applied for SOS and EOS, separately. First, the perpendicular distance of each threshold (first derivative, percentile 10%, percentile 20%, percentage 10% and percentage 20%) and smoothing methods (SavGol, HANTS, Polyfit) to the 1:1 line was calculated. The least mean distance of threshold and smoothing methods for each year was derived for SOS and EOS, separately. Through this analysis we obtained which threshold, or smoothing method, and which VI-based estimates of the date of SOS/EOS had the least difference with the date of SOS/EOS derived from GPP. For the second method, the mean difference between the SOS/EOS of GPP and the SOS/EOS of each VI for all threshold methods and smoothing functions were calculated for each year, separately. A positive difference illustrated an earlier SOS/EOS derived from VI timeseries, while a negative difference was indicative of a delay in SOS/EOS. Finally, a correlation analysis was performed to investigate the relationship between GPP and VIs. This analysis allowed us to understand if the coppicing in 2017 and drought in 2018 have an effect on the VI−GPP relationships and which VI can explain better this effect. We assumed the existence of a strong linear relationship between VIs and GPP as reported for different ecosystems [44,58]. Furthermore, this analysis was performed using original, not fitted, time-series data of both GPP and VIs to exclude the effect of filtering methods on the VI−GPP relationship. The goodness of the linear model was tested by the following statistics: the coefficient of determination (R 2 ) and root mean square error (RMSE).

GPP Temporal Variability
The SRC plantation showed very high plant photosynthesis, as indicated in Figure 2a by cumulative daily GPP values. The plant photosynthesis varied largely depending on environmental conditions and management activity (Figure 2). The seasonal and interannual variability of photosynthesis was associated mainly with soil water content (SWC, Figure 2b). The 2016 growing season, which corresponds to the third and the last year of the third rotation, showed the highest values of daily GPP. The constant and high values of SWC, for most of the year, did not limit the carbon photosynthesis uptake for that year. In contrast, the photosynthesis in the 2018 growing season, which had the hottest and driest summer, was reduced by a sharp decrease of SWC at the end of the summer and a low SWC until September. Compared to 2016 and 2018, an earlier recharge of SWC was observed after the summer in 2017 (Figure 2b). This increase of SWC arose GPP activity at the end of the 2017 growing season. Air temperature and solar radiation were similar between the three growing seasons (Figure 2a).
The coppicing strongly reduced the photosynthetic activity of the SRC plantation in spring 2017 (Figure 2a). The maximum GPP was only half that in other years (15-30 g C m −2 d −1 ).

Describing SOS and EOS Using VIs
The three filtering methods applied to the raw data time-series produced different seasonal patterns (Figure 3). For GPP, the filters HANTS and SavGol described well the GPP pattern at SOS and EOS, while Polyfit did not perform well at EOS (Figure 3a). Furthermore, the HANTS filter captured best the GPP amplitude, and the SavGol filter followed best the interannual variability in GPP (Figure 3a). However, the estimates of phenological metrics of GPP showed a very small variability between years: day of the year (doy) 115 (SD ± 6) for SOS and doy 313 (SD ± 8) for EOS.
The For VIs, the three filtering methods showed similar good performance in reproducing seasonal and interannual variability of GPP time-series at SOS (Table 3). The accuracy varied largely among VIs (Table 3, RMSE = 0.09-14.33 mol m −2 s −1 ). Better performance in predicting GPP pattern at the SOS than at the EOS (Table 3, Table 1). RMSE-root mean squared error, in mol m −2 s −1 . Scatter plots in Figure 4 shows the comparison between SOS derived from VIs and from GPP, grouped by year for threshold and filtering methods. The results of the same analysis applied to EOS are represented by scatter plots in Figure 5. A good performance in predicting SOS of GPP for all threshold methods and years was obtained for MTCI (Figure 4), points were the closest to the 1:1 line. For MCARI, NDVI, EVI and ChlRedEdge, the first derivative method failed in estimating SOS of GPP in 2016 for all filtering methods (Figure 4). Points for EOS for MTCI and EVI were the closest to the 1:1 line ( Figure 5) indicating good performance in predicting EOS of GPP. Except for the first derivative method in 2016 and 2018, MCARI showed good performances too. Among other VIs, NDVI, gNDVI, ChlRedEdge and SIPI showed a clear tendency to predict the EOS later than GPP.

RMSE
The analysis of the perpendicular distance to the 1:1 line confirmed that HANTS and Polyfit were the most suitable filtering methods for SOS, while for EOS it was not possible to identify only one method (Table 4). For SOS, the MTCI showed the minimum distance to 1:1 line (3.10 days (SD ± 3)), followed by ChlRedEdge (4.7 days), NDVI (5.08 days) and EVI (5.18 days). However, ChlRedEdge and NDVI showed very high variability (SD ± 4.23 and SD ± 4.25, respectively) among the different filtering methods (Figure 3 and Table 4). All other VIs showed larger distances, indicative of poor performance in predicting the SOS of GPP phenology. For EOS, MTCI and EVI showed a small distance to 1:1 line (10.41 days (SD ± 4.18) and 11.78 days (SD ± 5.48), respectively). Bigger distance values were found for other VIs (Table 4). Furthermore, the distance to the 1:1 line was VI-dependent and varied strongly with the filtering method (Table 4). However, the distance analysis confirmed that MTCI predicted well the SOS in 2016 (2.12 days (SD ± 1.81)) and the EOS in 2016 and 2018, the year with limiting SWC at the end of the growing season (Figure 3; 11.03 days (SD ± 2.49) and 6.59 days (SD ± 1.99), respectively). MCARI predicted well GPP at SOS in 2017 (the year of the copping; 1.37 days (SD ± 0.29), gNDVI (1.40 days (SD ± 0.84) at the SOS in 2018, while NDVI showed the best accuracy at the EOS in 2017 (7.68 days (SD ± 2.49; Figure 3, Table 4).
The 20% and 10% percentile were the best threshold for estimating GPP metrics for SOS. The distance between the EOS of GPP and EOS of VIs differed the least using a 10% percentile threshold (Table 5).   Bar plots in Figure 6 show the difference between the SOS/EOS derived from GPP and the SOS/EOS derived from VIs for each year averaged for filtering methods. A positive difference illustrated an earlier SOS/EOS derived from VI time-series, while a negative difference was indicative of a delay in SOS/EOS. SIPI predicted more than one month later the SOS of GPP and about two months earlier the EOS of GPP. All VIs (except SIPI and MCARI) estimated the SOS earlier than GPP in 2016. MTCI showed the smallest difference with the SOS of GPP. This tendency to predict earlier the SOS of GPP was also evident for all VIs, except for MTCI, MCARI and SIPI, in 2018. PSSR predicted about 2 days earlier the SOS of GPP. For the 2017 growing season, the year of the copping, all VIs, except ChlRedEdge and gNDVI, delayed the SOS of GPP. NDVI showed the smallest difference with the SOS of GPP. Clearly, all VIs (except SIPI and MCARI in 2016, and SIPI in 2017) predicted later the EOS than the GPP for all years. The EOS predicted by EVI and MCARI, showed a smaller difference with the EOS of GPP in 2017 and 2018. PSSR differed little to EOS of the GPP in 2016. The VIs performed differently in predicting daily GPP of the SRC plantation across the three growing seasons ( Figure 6 and Supplementary material S1). In addition, performances in predicting GPP also varied greatly across the years. MTCI performed best in predicting GPP for the 2016 and 2017 growing seasons ( Figure 6) and for all three growing seasons together (R 2 = 0.52, RMSE = 22.53 mol m −2 s −1 ). In contrast, ChlRedEdge performed slightly better for the 2018 dry growing season (R 2 = 0.59, RMSE = 22.86 mol m −2 s −1 ). EVI showed a slightly lower performance than ChlRedEdge (R 2 = 0.41, RMSE = 25.77 mol m −2 s −1 ) while NDVI showed a poor performance in predicting GPP (R 2 = 0.31, RMSE = 25.69 mol m −2 s −1 ). All other VIs had limitations in predicting GPP seasonality (see Supplementary material S1).

Discussion
This analysis explores for the first time the potential of using vegetation indices derived from the Sentinel-2 satellite for describing carbon flux phenology and productivity of SRC plantations.

GPP Phenology of SRC Plantations
The SRC plantation exhibited strong photosynthetic activity, in terms of GPP, until late in the season (Figure 2a). This was in line with previous studies on the same plantation reporting a net carbon uptake until the end of September [11,19], and on poplar trees [20]. Contrary to SRC plantations, a marked decline of net carbon uptake in early mid-September was observed in deciduous species [18]. The different growth strategy between SRC plantations and other deciduous species could be explained by a high allocation of nitrogen in leaves until late in the growing season, as reported for the same SRC plantation in . In addition, in this SRC plantation, the production of new leaves quickly decreases around September [100].
The year-to-year variations in water and temperature regimes (Figure 2b) affected the phenology of the plantation (Figures 2 and 3). Low temperatures in spring 2016 delayed the start of the growing season in that year. Low SWC at the end of season reduced GPP and anticipated the senescence in 2017 and 2018 ( Figure 2). However, the SWC recharge after the summer in 2017 contributed to a delay at the end of the season in 2017 compared to 2018. Poplars show a high photosynthetic performance when SWC is not a limiting factor and GPP is controlled by a decrease in solar radiation and air temperature [17]. A reduction of photosynthetic activity under dry conditions might be caused by a decrease in leaf production and stomatal closure during dry periods [17].
Coppicing in 2017 had a distinct impact on GPP seasonality and interannual variability of the plantation (Figure 2). A delay in the start of the season was observed after copping and was documented for the same plantation in the previous study of Vanbeveren et al. (2016). After coppicing, poplars invest energy in the production of new buds from stumps, a process that delays the production of news leaves and the start of photosynthesis. However, this delay in the season start was counterbalanced by rapid growth and resprouting in subsequent growing seasons [101].
We found a delay of EOS in 2017 compared to 2018 and 2016, but there was no evidence that this delay was directly linked to copping. Based on this, and on previous studies reporting that coppicing does not affect autumn senescence [21], it is possible that SWC limitations might rather have contributed to the delay of EOS observed in 2017.

An Optimal Spectral Proxy for Describing GPP Phenology of SRC Plantations
Overall, our results revealed that MTCI was the most suitable proxy for describing the GPP phenology (i.e., it exhibited the lower distance in identifying the SOS-and EOS-dates (Tables 4 and 5 Table 4). The estimates of phenological metrics derived from MTCI showed small differences with GPP phenology for all years anticipating the SOS and delaying the EOS. Among all other tested VIs SIPI and PSSR failed to reproduce the seasonal GPP pattern and to represent SOS and EOS of the GPP. We found that also EVI provided quite good estimates of SOS and EOS and seasonality of GPP, albeit in all cases lower performance than MTCI. NDVI, ChlRedEdge and MCARI performed a bit weaker than EVI, at both SOS and EOS. NDVI performed weaker than EVI, which performed weaker than ChlRedEdge in simulating GPP seasonal variability (Figure 7).
MTCI is well known to be a reliable proxy of canopy chlorophyll content [35,52]. It has been successfully used as a proxy of GPP seasonal variability in several biomes [33,53,[102][103][104][105]. Changes in canopy chlorophyll content play an important role in controlling ecosystem photosynthetic activity, i.e., GPP [49,105]. We observed an exceptionally good representation of the seasonal GPP pattern by MTCI ( Figure 3, Table 3, Figure 7) [34,35,106]. This tight relationship between GPP and MTCI has been attributed to the fact that MTCI does not saturate at higher chlorophyll contents [69,107].
As expected, NDVI was not the best proxy for GPP (Tables 4 and 5; Figure 7). NDVI index-such as the widely used "greenness indices"-has been demonstrated to be a good proxy of the fraction of absorbed photosynthetically active radiation (i.e., fAPAR) and LAI [108,109]. Changes in LAI among phenophases affect the shape of the relationship between NDVI and LAI that has almost linear trends during the growing season and declining during the senescence period [110,111]. NDVI shows limitations in describing seasonal LAI variations of dense canopies because it responds mainly to variations in the red and it is relatively insensitive to variations in the NIR [37,112]. These SRC plantations reached an LAI of about 10 in 2016 (data not shown), the last year of the third rotation, when the canopies were fully developed. The very high LAI in 2016 can explain the poor performance of NDVI for the SRC plantations. LAI of natural mature forests varies between 2 and 7 [37] and it could explain the better performance of NDVI in describing the phenology of deciduous forests. The same considerations can explain also the poor performance of other VIs calculated using the NIR reflectance (Table 1: SIPI, PSSR and gNDVI). EVI is assumed to relate to green cover and, therefore, be a useful proxy for canopy phenology [113][114][115][116]. Other studies also reported that EVI well represents the seasonal patterns of GPP [42,117]. EVI was formulated to be less sensitive to atmospheric effects and reflection from soil than other greenness indices such as NDVI, by incorporating blue spectral wavelengths [118] to differentiate soil from vegetation [119].
The chlorophyll red-edge index (ChlRedEdge) is well recognized as a good proxy of chlorophyll content [120] and has been demonstrated to be a valuable descriptor of the seasonal changes in GPP [121]. Also in our study, ChlRedEdge showed good accuracy in predicting SOS and seasonality of GPP, albeit less than MTCI (refer to Tables 4 and 5; Figure 7). Recently, it has been shown that the ChlRedEdge index can be a reliable proxy for leaf senescence as well for deciduous trees (i.e., European beech, pedunculate oak and silver birch) in the temperate zone [122]. However, this was not the case in our study, as ChlRedEdge did not provide a good estimation of EOS. It is important to note, however, that Mariën et al. (2019) focused on the onset of autumn senescence as identified by the breakpoint technique, whereas we used low percentiles and percentages, which correspond with the later stages of senescence, which could explain the different outcome between both studies. Coppicing strongly affects the canopy structure (i.e., leaf angle and distribution, orientation, leaf density) of the plantation modifying leaf light interception. These changes in canopy structure can strongly modify the spectral response of SRC plantations compared to natural deciduous forests. In addition, canopy structure can have a strong impact on the estimation of canopy nitrogen content [123]. During senescence, SRC plantations allocate a larger amount of nitrogen in leaves than deciduous forests. It could explain the different spectral responses of the SRC plantation compared to deciduous forests.
Sentinel-2 data are now systematically generated and the Level-2A product provides top of canopy reflectance images globally (https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/producttypes/level-2a). This permits to test the use of Sentinel-2 vegetation indices, in particular the MTCI, over other SRC plantations and different biomes. The recent ESA Sentinel-3 mission offers a unique set of data at the global scale with an improved spectral resolution of about 10 nm that may help to better understand plant traits and plant physiological status. The Sentinel-3 OLCI sensor covers the chlorophyll absorption region between 600 nm and 677 nm and the red-edge region between 697 nm and 755 nm. Sentinel-3 Level-1 OLCI products provide the so-called OLCI Terrestrial Chlorophyll Index (OTCI) which is a chlorophyll index that has the same formula of MTCI but it is calculated using the OLCI bands. OTCI is available at a spatial resolution of 300 m instead of the 20 m of Sentinel-2. Thus, the use of an approach similar to the one based on MTCI and presented here could contribute to the understanding of SRC phenology and GPP at a global scale.

Quality of Satellite Raw Data, Filtering and Threshold Methods
The filtering methods HANTS and Polyfit reproduced a better method for the seasonal and interannual variability of GPP of the SRC plantation ( Figure 3; Table 3) than the SavGol method. Furthermore, both filtering methods were suitable for identifying the start and the end of carbon flux uptake (Table 4), as reported in previous studies [124]. Atkinson et al. (2012) [125] reported that fitting methods where each parameter of the fitting function is estimated by the least-square method are more robust. However, there is not a unique method that always performs better than another for filtering VI time-series [64].
The differences between the different filtering methods in reconstructing VI time-series could be related to their ability in detecting clouds and cloud shadow. The presence of clouds changes the optical properties of the signal. The cloud shadow weakens the light intensity and thus the reflectance of all bands is underestimated [64]. Furthermore, while the presence of clouds causes an underestimation of VIs, cloud shadow can lead to an overestimation of NDVI [64]. A limited description of the EOS in all years by most of the VIs could be also related to the large presence of clouds at the end of the summer and in autumn when senescence started ( Figure 3, Table 4).
Another point worth mentioning relates to the accuracy of cloud masks used in the screening of clouds and cloud shadows. In this study, we used postprocessed Sentinel-2 data by Sen2Cor algorithm. The cloud classification generated by Sen2Cor shows limits in accurately identifying the cloud shadows. This classification does not reduce the noise of the time-series data but it produces a misleading classification of the satellite data quality. Therefore the remaining data with low quality but classified as good pixels might affect the patterns of fitted VI time-series ( Figure 3, Table 3). The number of high-quality data of the timeseries and their distribution in time affects the accuracy of the fitting functions [63]. A comparison of filtering methods is advisable, given their contrasting performance in the literature [64,89,125].
In addition to the quality of VI time-series and the robustness of the filtering functions, the accuracy of the phenophases also depends on the chosen detection method [126], i.e., the selected dates of SOS-and EOS occurrence depend on the thresholds that are chosen to represent them (Figures 4 and 5). In this study, we found that for SOS, the dates at which the 10% and 20% percentile occurred showed very good agreement between VIs and the GPP data, while for EOS the dates of the 10% percentile came out as the best threshold ( Table 5). The percentile method focuses on the position of the proxy for a fixed percentage (e.g., 10% and 20%). Several previous studies also obtained good results using such low percentiles as thresholds for determining the dates of SOS and EOS [127][128][129][130][131]. Furthermore, a threshold around 20% is good for estimating the start of greening [132]. Our analysis confirmed also that the extraction of phenological metrics from VI time-series is highly uncertain [65][66][67]. Combining multiple methods is preferred especially for estimating EOS dynamics ( Table 5). The first derivative method is very sensitive to the noise in the signal and the temporal filtering and it cannot represent short growing seasons well, especially when the increase and decrease in the annual time-series occur rapidly and abruptly [133]. The percentage method is sensitive to the minimum and maximum values that may be affected by noise in the signal.

Conclusions
This study confirmed that VIs can be suitable estimators of the GPP phenology of SRC plantations and are able to detect the effect of coppicing on GPP. Among all investigated vegetation indices, MTCI showed the best fit with the GPP phenology. Additionally, the choice of the filtering method and the thresholds