1 Land surface phenology from VEGETATION and 2 PROBA-V data. Assessment over deciduous forests 3 4 Kevin Bórnez a, b, *, Adrià Descals a, b Aleixandre Verger a, b, Josep Peñuelas a, b 5 a CREAF, Cerdanyola del Vallès 08193, Catalonia, Spain 6 b CSIC, Global Ecology Unit, Cerdanyola del Vallès 08193, Catalonia, Spain 7 * Corresponding author. E-mail address: k.bornez@creaf.uab.cat (k. Bórnez) 8 ABSTRACT 9 Land surface phenology has been widely retrieved although no consensus exists on the optimal satellite 10 dataset and the method to extract phenology metrics. This study is the first comprehensive comparison of 11 vegetation variables and methods to retrieve land surface phenology for 1999-2017 time series of Copernicus 12 Global Land products derived from SPOT-VEGETATION and PROBA-V data. We investigated the sensitivity 13 of phenology to (I) the input vegetation variable: normalized difference vegetation index (NDVI), leaf area index 14 (LAI), fraction of absorbed photosynthetically active radiation (FAPAR), and fraction of vegetation cover 15 (FCOVER); (II) the smoothing and gap filling method for deriving seasonal trajectories; and (III) the method to 16 extract phenological metrics: thresholds based on a percentile of the annual amplitude of the vegetation variable, 17 autoregressive moving averages, logistic function fitting, and first derivative methods. We validated the derived 18 satellite phenological metrics (start of the season (SoS) and end of the season (EoS)) using available ground 19 observations of Betula pendula, B. alleghaniensis, Acer rubrum, Fagus grandifolia, and Quercus rubra in Europe 20 (Pan-European PEP725 network) and the USA (National Phenology Network, USA-NPN). The threshold-based 21 method applied to the smoothed and gap-filled LAI V2 time series agreed best with the ground phenology, with 22 root mean square errors of ~10 d and ~25 d for the timing of SoS and EoS respectively. This research is expected 23 to contribute for the operational retrieval of land surface phenology within the Copernicus Global Land Service. 24 Keywords 25 Land surface phenology, SPOT-VEGETATION, PROBA-V, leaf area index, ground observations 26 1. Introduction 27 Phenology is the study of the timing of recurrent biological and seasonal events and their 28 biotic and abiotic factors (Beaumont et al., 2015). Studies of plant phenology focus on how 29 these events and factors are influenced by seasonal and interannual variations in climate and 30 how they modulate abundance and diversity (Beaumont et al., 2015). Phenology is, moreover, 31 key to control physicochemical and biological processes, especially albedo, surface 32 roughness, canopy conductance and fluxes of carbon, water and energy (Peñuelas et al., 2009; 1 33 Richardson et al., 2013). Phenological metrics are thus relevant parameters for modeling land 34 surface processes and the global carbon cycle (Wu et al., 2014). 35 Phenological metrics are estimated based on ground observations and data derived from 36 satellites. Ground observations provide accurate timing of vegetation phenophases but cannot 37 cover continuously large-scale areas (Garrity et al., 2011; Yu et al., 2017). Satellite sensors 38 with moderate spatial resolutions, including AVHRR, MODIS, MERIS, SPOT- 39 VEGETATION and PROBA-V, provide long-term time series of daily observations that allow 40 improving the characterization of land surface phenology on a global scale (Verger et al., 41 2016, Atkinson et al., 2012; Zhang et al., 2004). However, the noise in the data and missing 42 observations mainly due to cloud contamination may induce significant uncertainties in the 43 estimation of phenological metrics (Kandasamy et al., 2013; Verger et al., 2013). The 44 literature shows a broad variety of time-series processing methods designed to reconstruct 45 gap-filled vegetation seasonal trajectories from noisy satellites signals. This includes the best 46 index slope method (Viovy et al., 1992), mean filters (Reed et al., 1994), moving-window 47 filters (Sweets et al. 1999), asymmetric Gaussian functions (Jönsson & Eklundh, 2002), 48 Savitzky–Golay filters (Chen et al., 2004) or the Whittaker smoother (Eilers, 2003). However, 49 no single method always performs better than others for smoothing vegetation time series (Cai 50 et al., 2017) and their performance vary spatially and temporally with land surface conditions 51 and cloud influence (Atkinson, et al. 2012; Kandasamy et al., 2015). 52 A broad variety of statistical methods have been designed to extract phenological metrics 53 from satellite time series. Metrics typically include the start of the season (SoS), the end of 54 the season (EoS), the timing of maximum growth and the length of the growing season (LoS) 55 (Reed et al., 1994; Zhang et al., 2004). de Beurs and Henebry (2010) provided a 56 comprehensive review of the exiting phenology retrieval approaches that can be classified in 57 four main categories: thresholds and percentile based methods (Atzberger et al., 2011; Verger 58 et al., 2016), moving averages (Reed et al., 1994), first derivatives (White et al., 2009) and 59 fitted models (de Beurs and Henebry, 2005). White, et al. (2009) compared ten different 60 phenology retrieval methods applied to AVHRR NDVI in North America and found large 61 discrepancies of up to two months in the detection of the SoS. 62 In addition to the sensitivity to the smoothing and phenological extraction algorithm, the 63 derived phenological metrics are also dependent on the sensor, spatial and temporal resolution, 64 processing chain, and satellite data set. The satellite-derived spectral vegetation indices (e.g. 2 65 the Normalized Difference Vegetation Index (NDVI)) vary in their strength of phenological 66 prediction across sites and plant functional types (Wu et al., 2014). Unlike previous studies 67 based on vegetation indices, the present study aimed to characterize the phenology not only 68 with NDVI but also with biophysical variables: the leaf area index (LAI), the fraction of 69 absorbed photosynthetically active radiation (FAPAR), the fraction of vegetation cover 70 (FCOVER). We used NDVI version V2.1 (Toté et al., 2017), LAI, FAPAR and FCOVER V1 71 (Baret et al., 2013) and V2 (Verger et al., 2014) time series derived within the Copernicus 72 Global Land Service (CGLS) from SPOT-VEGETATION and PROBA-V data. Verger et al. 73 (2017) showed that the phenology derived from the interannual climatology of LAI V1 74 improved other existing products including MODIS-EVI when compared to ground 75 observations for the average date of the SoS and EoS. However, their study was limited to the 76 baseline LAI phenology as derived from a single extraction method. This paper is a 77 continuation of the previous paper by Verger et al. (2017) and we address now the interannual 78 variation of the yearly phenology, the impact of the input vegetation variable and the 79 phenological extraction method. Further, we incorporate LAI, FAPAR and FCOVER V2 that 80 improved continuity (no missing data in V2) and smoothness as compared to V1. 81 Our study had two main objectives: to select the best biophysical variable or vegetation 82 index for estimating phenological metrics on a global scale within the portfolio of the CGLS 83 vegetation products (NDVI, LAI, FCOVER or FAPAR) and to define the method that best 84 matched the ground data. 85 2. Materials and methods 86 2.1. Phenological ground observations 87 Ground-based phenological data from PEP725 and USA-NPN were examined, focusing 88 on the dates of leaf out and leaf senescence for Betula (birch) in Europe (Figure 1a) and the 89 USA and for Quercus (oak), Fagus (beech) and Acer (maple) in the USA (Figure 1b). These 90 genera were chosen because they are present in both Europe and the USA and have large 91 numbers of records in the combined data set. 92 The PEP725 Pan-European Phenology database (Templ et al., 2018) (www.pep725.eu) has 93 complete records from 1990 to the present. The phenophases defined in PEP725 are based on 94 BBCH (Biologische Bundesanstalt, Bundessortenamt and Chemical industry) code (Meier, 95 2009). We used the phenophases corresponding to the first visible leaves (BBCH 11) as the 3 96 reference for the timing of SoS and the date corresponding to 50% of leaves with autumn 97 coloration (BBCH 94) for the timing of EoS. 98 The USA National Phenology Network was established in 2007 to collect, store and share 99 historical and contemporary phenological data on a North American scale (Schwartz et al., 100 2012). The data are freely available at https://www.usanpn.org/. This network provides 101 measurements of several phenophases. We used the phenophases “leaves” which corresponds 102 to first visible leaves and “increasing leaf size” for SoS and “colored leaves” for EoS. 103 We discarded ground sites with less than four yearly measurements to obtain consistent 104 data records over the time series. We used ground-site located pixels. The spatial 105 heterogeneity hamper comparing ground-based phenology for individual plants with satellite 106 phenology at a resolution of 1 km. We filtered the ground sites located in agricultural or urban 107 areas using high-resolution images from Google Earth (https://earth.google.com/) and the 108 ESA Land Cover Map (CCI-LC) (http://maps.elie.ucl.ac.be/CCI/viewer/index.php). In 109 Europe, we used only the points with forest coverages >5% using a tree cover map for 110 European forests (Brus et al., 2012). 111 112 Fig. 1. Location of selected phenological ground observations at (a) PEP725 sites for Silver birch (Betula 113 pendula) and (b) USA-NPN sites for Red maple (Acer rubrum), Yellow birch (Betula alleghaniensis), American 114 beech (Fagus grandifolia) and Northern red oak (Quercus rubra). 115 2.2. Satellite time series 116 The time series of satellite imagery used for estimating the phenological metrics were from 117 the SPOT-VEGETATION (1999-2013) and PROBA-V satellites (2014-2017) with spatial 4 118 resolutions of 1 km and temporal frequencies of 10 d. In particular, we used LAI, FAPAR, 119 FCOVER V1 (Baret et al., 2013) and V2 (Verger et al., 2014) and NDVI V2.1 (Toté et al., 120 2017) products generated within the CGLS 121 (https://land.copernicus.eu/global/themes/vegetation) (Table1). 122 Time series of LAI, FAPAR and FCOVER V1 and NDVI V2.1 contained frequent 123 disturbances caused by residual cloud contamination, atmospheric variability, snow and bi- 124 directional effects (Figure 2). These time series required the application of smoothing and gap 125 filling techniques to generate consistent and gap filled seasonal trajectories (section 2.3) 126 before the extraction of phenological metrics. The LAI and FCOVER V2 products improved 127 over the V1 products in terms of temporal consistency and continuity (Figure 2). The V2 128 algorithm included multi-step data filtering, smoothing and gap-filling techniques that 129 rendered the products suitable for phenological estimation without additional pre-processing. 130 Filtering of outliers is based on an upper envelope approach, the pixel climatology 131 (interannual mean) is used to fill missing data and a Savitzky-Golay filter is used for the 132 smoothing (Verger et al., 2014). 133 Table 1. Algorithm principles of NDVI V2.1, and LAI, FAPAR, FCOVER V1 and V2. NDVI LAI, FAPAR, FCOVER LAI, FAPAR, FCOVER Version 2.1 Version 1 Version 2 Inputs Top of the canopy Nadir normalized TOC TOC reflectances in the red, (TOC) reflectances in reflectances in the red, NIR NIR and SWIR spectral the red and near and short-wave infrared bands, and cosine of the 3 infrared (NIR) (SWIR) spectral bands, and angles of sun and view spectral bands cosine of the sun zenith directions angle at 10:00 local time Temporal 10 d compositing 30 d compositing period Adaptive compositing composition period. The maximum with Gaussian weighting within 15 and 60d semi- NDVI value in the (minimum of two valid periods defined by the composition window observations) availability of 6 valid is retained. Nominal dates: 3rd, 13th observations at each side of Starting date of and 21-24th day of each the date being processed composition:1st, 11th month Nominal dates: 10th, 20th and 21th day of the and last day of the month month Temporal Not applied Not applied Multi-step filtering, smoothing and temporal smoothing and gap filling gap-filling 5 134 135 Fig. 2. LAI (V1 and V2) and NDVI (V2.1) time series for the PEP725 site 5449 (50º42'20.49"N, 13º46'59.55"E) 136 representative of birch forest in Europe. 137 2.3. Smoothing methods 138 We tested several smoothing methods for reducing noise and reconstructing gap filled 139 seasonal trajectories from CGLS time series (Eerens and Haesen, 2015): 140  WHITTAKER smoother (Atzberger et al., 2011): It minimizes a cost function 141 describing the balance between fidelity (quadratic difference between estimates and 142 actual observations) and roughness (quadratic difference between successive 143 estimates). 144  BISE (Best Index Slope Extraction) (Viovy et al., 1992): It retains the good 145 observations in a local window and replaces missing or eliminated suspect values by 146 linear interpolation. 147  MEAN: A linear interpolation is first applied to fill missing data. A running mean filter 148 with a sliding window of 50 d length is then applied. 149  SWETS method (Swets et al., 1999): A linear interpolation is first applied to fill missing 150 data. A weighted linear regression over a local window is then applied. 151 2.4. Methods for extracting phenological metrics 152 We tested four state of the art methods to extract phenological metrics from CGLS time 153 series (Figure 3): 154  Thresholds based on a pixel percentile value: SoS is defined as the day of the year 155 (DoY) when a vegetation variable exceeds a particular threshold. EoS is defined as the 156 DoY when an index remains below a particular threshold. We established dynamic 6 157 thresholds per pixel based on a percentile of the annual amplitude of the vegetation 158 variable (Verger et al., 2016). The selected percentiles were determined based on the 159 comparison with available ground measurements. We tested the 20th, 30th, 40th and 160 50th percentiles of the annual amplitude for SoS, and the 30th, 40th, 50th and 60th 161 percentiles for EoS. This method is also the basis of SPIRITS phenological approach. 162  Autoregressive moving average: A moving average is first computed at a randomly 163 chosen time lag (Ivits et al., 2009). We tested time lags from 50 to 150 d and selected 164 a time lag of 100 d based on the comparison with ground measurements. SoS and EoS 165 are then defined as the DOY when the moving average curves cross the original curve 166 of the vegetation variable. 167  First derivative: SoS is defined as the DoY of the maximum increase (maximum first 168 derivative) in the curve (Tateishi & Ebata, 2004). EoS is defined as the DoY of the 169 maximum decrease in the curve. 170  Logistic function: SoS is defined as the DoY of the first local maximum rate of change 171 in the curvature of a logistic function fitted to the time series (Zhang et al., 2003). EoS 172 is defined as the DoY of the first local minimum rate of change in the curvature. 173 Fig. 3. Schematic representation of SoS (on the left of the peak) and EoS (on the right of the peak) retrieved with 174 the four methods for the PEP725 site 4959 (50º42'20.49"N, 13º46'59.55"E) for 2011. The black circles 175 correspond to the original LAI data at a 10-d frequency, and the green line corresponds to the data interpolated 176 at daily steps, which is used for phenological estimation. 177 178 7 179 2.5. Methodological approach 180 The several satellite-derived vegetation variables, smoothing methods and phenological 181 extraction approaches lead to a large number of combinations. We sequentially investigated 182 the impact of smoothing, variable and extraction method based on an initial set of modalities. 183 The initial modalities were defined a posteriori based on the analysis of all the combinations: 184 1. Sensitivity analysis of the smoothing method: We used LAI V1 as input dataset and 185 the percentile phenology method. 186 2. Sensitivity analysis of the vegetation variable: We used the SWETS smoothing 187 method and the percentile phenology method. 188 3. Sensitivity analysis of the method to extract phenological metrics: We used the LAI 189 V2. Note that in this case the application of a smoothing method is not required 190 because LAI V2 is already smoothed and gap-filled. 191 The analysis 1 and 2 were carried out in Europe and for the validation we used ground 192 measurements of Betula Pendula, which showed a greater latitudinal distribution. The 193 analysis 3 was performed at the global scale. For the sensitivity analysis 1 and 2, we used the 194 Software for the Processing and Interpretation of Remotely Sensed Image Time Series 195 (SPIRITS) (Eerens et al., 2014; Eerens and Haesen, 2015). For the sensitivity analysis 3, we 196 used Google Earth Engine (GEE) (https://earthengine.google.org) which allowed 197 implementing dedicated algorithms while only the threshold method is available in SPIRITS. 198 The input 10 d time series were linearly interpolated at daily steps before phenological 199 retrieval. For SPIRITS, the precision of phenological estimates is limited by the frequency of 200 the input time series (10 d in our case) (non interpolation). For pixels with multiple growing 201 seasons, we computed the phenological metrics for the growing season having the highest LAI 202 amplitude. 203 The agreement between metric estimates from satellite imagery and ground-based 204 measurements was quantified using the slope of the linear regression, the Pearson correlation 205 coefficient (R), bias, i.e. the average difference between the satellite-derived phenology and 206 the observed date (a positive bias indicated that SoS and EoS occurred later than the observed 207 leaf out and autumnal coloring, respectively), the absolute bias and root mean square error 208 (RMSE) calculated (e.g. for SoS) as: 8 (1) 209 where n is the number of samples 210 3. Results 211 3.1. Sensitivity analysis of the smoothing method 212 The SWETS method performed the best for the reconstruction of seasonal trajectories and 213 the estimation of phenological metrics:16 d in terms of RMSE and 1 d in terms of bias for the 214 timing of SoS, and 34 d (RMSE) and 12 d (bias) for the timing of EoS (Table 2). 215 Table 2. Statistics of comparisons between the derived phenology metrics from reconstructed LAI V1 time series 216 using the different smoothing methods and PEP725 ground measurements (Betula pendula) for the SoS and EoS. 217 The SPIRITS percentile method was used to extract the phenology using 30% of LAI amplitude for the SoS and 218 40% for the EoS. * Significant correlations with p<0.05 (**p<0.001). Metric Method RMSE (d) Bias (d) Abs. Bias (d) R Slope SoS WHITTAKER 23.80 -8.16 16.76 0.41 0.76 "visible BISE 27.22 -7.52 19.84 0.43* 0.89 leaves" MEAN 23.98 -5.84 18.43 0.39 0.74 (n=359) SWETS 16.44 -1.38 17.39 0.49** 0.67 EoS WHITTAKER 44.56 13.51 31.65 -0.19 -0.12 "colored BISE 49.45 19.50 37.25 -0.03 -0.10 leaves" MEAN 44.10 13.18 30.54 -0.08 -0.01 (n=359) SWETS 34.28 12.01 28.96 -0.15 -0.19 219 220 3.2. Sensitivity analysis of the vegetation variable 221 The best agreement with ground measurements for the timing of SoS was found for the 222 30% threshold of LAI amplitude, 40% of FCOVER, 50% of FAPAR and 50% of NDVI. The 223 best metric definitions for the EOS were based on 40% of LAI and FCOVER amplitudes and 224 50% for FAPAR and NDVI. Phenological metrics derived from V2 time series improved over 225 V1 for all variables (compare V1 and V2 statistics in Table 3). 226 227 9 228 Table 3. Statistics of comparisons between the derived phenologies from LAI, FCOVER, FAPAR V1 and V2 229 and NDVI V2.1 time series and PEP725 ground measurements (Betula pendula) for the start of the season (SoS) 230 and end of the season (EoS). The SPIRITS percentile method was used to extract the phenology with specific 231 thresholds per vegetation variable and phenological metric. *mark indicates significant correlations with p<0.05 232 (** indicates p<0.001). Metric Index Version Definition RMSE (d) BIAS (d) Abs. BIAS (d) R Slope SoS LAI V1 30% 16.44 -1.38 17.39 0.49** 0.67 (n=359) V2 30% 12.49 1.65 10.22 0.62** 0.78 FCOVER V1 40% 17.25 -4.44 23.54 0.52** 0.74 V2 40% 13.95 -4.46 11.21 0.54** 0.95 FAPAR V1 50% 31.40 -18.35 31.23 0.48** 0.85 V2 50% 23.01 -14.26 13.94 0.57** 0.88 NDVI V2.1 50% 20.70 -13.18 15.38 0.25 0.82 EoS LAI V1 40% 34.28 12.01 28.96 -0.15 -0.19 (n=359) V2 40% 32.72 -6.18 21.12 0.05 0.11 FCOVER V1 40% 30.69 -14.33 42.25 -0.13 -0.20 V2 40% 25.35 6.42 20.34 0.26 0.38 FAPAR V1 50% 44.06 18.75 45.87 -0.23 -0.65 V2 50% 39.99 18.61 37.97 -0.04 -0.51 NDVI V2.1 50% 48.35 14.75 34.58 0.02 -0.25 233 234 The best performances for SoS were obtained using the LAI and FCOVER V2 time series 235 (Figure 4a), with RMSEs of ~12 and 14 d, respectively (Table 3). In contrast, RMSEs were 236 ~21 and 23 d for NDVI V2.1 and FAPAR V2, respectively (Table 3), and contained many 237 outliers (Figure 4a). SoS was slightly underestimated for all cases, except when using LAI 238 (Figure 4a). EoS had a higher RMSE (25-48 d) and a lower R (<0.3) (Table 3). The estimates 239 of EoS from the LAI and FCOVER V2 time series also agreed best with ground data although 240 no significant correlations were found (Table 3, Figure 4b). 10 241 242 Fig. 4. Boxplots of the bias errors for (a) SoS and (b) EoS estimated from the LAI, FCOVER, FAPAR V1 and 243 V2 and NDVI V2.1 time series minus the PEP725 ground measurements. An elongated boxplot indicates a 244 greater dispersion of the average bias. 245 3.3. Sensitivity analysis of the method to extract phenological metrics 246 The 30th percentile and SPIRITS applied to LAI V2 provided the best performances among 247 the different analyzed methods when compared both with USA-NPN and Europe-PEP725 248 measurements of the SoS (Figure 5, Table 4). The results for the 30th percentile and SPIRITS 249 were similar because both methods use the same definition of SoS based on the 30% threshold 250 of annual amplitude but the 30th percentile method slightly improved SPIRITS in terms of 251 precision (c.f. Figure 6a-6b, 7a-7b) and accuracy (RMSE of 9 vs 14 d for USA-NPN “leaves”, 11 252 9 vs 10 d for USA-NPN “increasing leaf size” and 11 vs 12 d for PEP725 “first visible leaves” 253 (Table 4)) because of the daily interpolation applied to the input 10 d data. The logistic 254 function also performed well (RMSE from 13 to 18 d) but provided slightly advanced SoS as 255 compared to ground measurements (bias from -5 to -9 d). The derivative method provided 256 showed higher RMSE (up to 28 d) and positive bias (up to 19 d). The moving-average method 257 performance the worst (RMSE from 30 to 59 d) and systematically advanced the timing of 258 SoS as compared to ground measurements (bias from -24 to -51d). 259 Table 4. Statistics of comparisons between LAI V2 derived phenology for the start of the season (SoS) and end 260 of the season (EoS), and the ground measurements for the various methods (percentiles, logistic function, 261 derivative and moving average). *mark indicates significant correlations with p<0.05 (** indicates p<0.001). Validation source Metric Method RMSE (d) Bias (d) Abs. Bias (d) R Slope USA-NPN SoS 30th percentile 9.19 1.59 6.77 0.81** 0.73 (n=462) "Leaves" SPIRITS 30% 13.75 0.58 7.05 0.60** 0.58 Logistic 14.14 -5.25 8.66 0.63** 0.60 Derivative 16.82 9.63 12.01 0.68** 0.82 Moving average 30.06 -24.33 23.95 0.59* 0.82 USA-NPN SoS 30th percentile 8.87 -2.20 6.11 0.83** 0.67 (n=158) "Increaising SPIRITS 30% 10.40 -3.46 6.18 0.78** 0.67 leaf size" Logistic 12.81 -7.24 6.86 0.74** 0.63 Derivative 10.04 4.22 7.73 0.83** 0.84 Moving average 34.49 -28.07 27.86 0.49* 0.70 PEP725 SoS 30th percentile 11.50 1.69 9.84 0.60** 0.90 (n=359) "visible SPIRITS 30% 12.49 1.65 10.22 0.62** 0.78 leaves" Logistic 17.96 -9.22 14.21 0.53* 0.71 Derivative 28.31 19.19 20.42 0.49* 0.88 Moving average 56.24 -51.35 54.27 0.50* 1.14 USA-NPN EoS 40th percentile 25.61 6.39 17.60 0.14 0.10 (n=241) "Colored SPIRITS 40% 30.79 -3.69 20.17 0.06 0.07 leaves" Logistic 30.70 21.80 23.66 0.10 0.15 Derivative 40.40 -10.11 26.70 0.02 0.01 Moving average 58.75 50.35 44.26 0.13 0.01 PEP725 EoS 40th percentile 27.69 -5.15 18.89 0.11 0.15 (n=359) "Colored SPIRITS 40% 32.72 -6.18 21.12 0.05 0.11 leaves" Logistic 30.15 11.91 24.29 0.03 0.28 Derivative 64.93 -43.96 48.49 0.00 -0.11 Moving average 54.95 37.81 52.55 0.01 0.21 262 263 12 264 Fig. 5. Boxplots of the bias error for the SoS estimated from LAI V2 minus the ground measurements at the 265 USA-NPN (a) and PEP725 (b) sites. An elongated boxplot indicates a greater dispersion of the average bias in 266 each method. 13 267 Fig. 6. Scatterplots between the SoS predicted from LAI V2 by the percentile method (a), SPIRITS (b), logistic 268 function (c), derivative (d) and moving average (e) compared with the ground phenology (USA-NPN “leaves”). 269 Values are given in DoY. Statistics of the comparison are indicated in Table 4. 270 Fig. 7. Scatterplots between SoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function (c), 271 derivative (d) and moving average (e) methods and ground phenology (PEP725 “first visible leaves”). Values 272 are given in DoY. Statistics of the comparison are indicated in Table 4. 14 273 The different methods provided poorer performances for the EoS (Table 4, Figure 8). The 274 best agreement with ground measurements was found for the 40th percentile method (RMSE 275 of 25 d for USA-NPN and 28 d for PEP725 “colored leaves” measurements). The logistic 276 function provided similar performances in terms of RMSE but slightly overestimated ground 277 measurements (bias from 12 to 22 d). The derivative method showed higher scattering (RMSE 278 from 40 to 65 d) and lower correlation than other methods (Table 4, Figure 9 and Figure 10). 279 The moving average retrievals showed a positive delay as compared to ground data (bias of 280 38 d for PEP725 and 50 d for USA-NPN). 281 282 Fig. 8. Boxplots of the bias error for the EoS estimated from LAI V2 minus the ground measurements at the 283 USA-NPN (a) and PEP725 (b) sites. An elongated boxplot indicates a greater dispersion of the average bias in 284 each method. 285 15 286 Fig. 9. Scatterplots between the EoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function 287 (c), derivative (d) and moving average (e) methods and ground phenology (USA-NPN “colored leaves”) for Acer 288 rubrum (green), Betula alleghaniensis (orange), Fagus grandifolia (blue) and Quercus rubra (yellow). Values 289 are given in DoY. Statistics of the comparison are indicated in Table 4. 290 291 Fig. 10. Scatterplots between the EoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function 292 (c), derivative (d) and moving average (e) methods and ground phenology (PEP725 “colored leaves”). Brown 293 for heterogeneous sites (forest cover <50%) and yellow for homogeneous sites (forest cover >50%). Values are 294 given in DoY. Statistics of the comparison are indicated in Table 4. 16 295 3.4. Spatial patterns of land surface phenology 296 Figure 11 shows the average timing of SoS, EoS and LoS phenophases at the global scale 297 by using the percentile method (30th and 40th percentiles for SoS and EoS, respectively) and 298 time series of V2 LAI (1999-2017). The LoS is estimated as the length of time between the 299 SoS and the EoS. The derived maps show consistent spatial patterns of the seasonality of 300 vegetation at the global scale which is driven by the distribution of latitudinal climatic 301 patterns, type of vegetation and topographic elements, among other factors (Verger et al., 302 2015; Zhang et al., 2004). The timing of SoS (Figure 11a) and EoS (Figure 11b) reflected a 303 broad variation in the range of values and their spatial pattern at the middle and high latitudes 304 of the Northern Hemisphere strongly depend on the thermal and photoperiod latitudinal 305 gradient (Verger et al., 2016). 306 The LoS (Figure 11c) in some ecoclimatic and biogeographic regions of transition such as 307 Sahel shows a broad range of variation from 10 to 200 d following the positive north-south 308 gradient of rainfall (Verger et al., 2016). On the contrary, in northern latitudes >50º the LoS 309 showed a limited range of variation from 10 to 60 d with shorter days at higher latitudes 310 following the latitudinal gradient of temperature and radiation (Verger et al., 2016). 17 311 312 Fig. 11. Global Map of average SoS (a), EoS (b) and LoS (c) derived from the LAI V2 time series (1999-2017) 313 and the threshold-based method. The continental areas in white are deserts and evergreen forests with limited 314 seasonality where phenology was not computed. 315 4. Discussion and conclusion 316 Land surface phenology (LSP) provides a synoptic view of vegetation dynamics and it can 317 substantially improve our macroecological knowledge and the representation of phenology in 318 earth-system models. Unlike previous studies limited to NDVI, we used three additional 319 biophysical variables: LAI, FAPAR and FCOVER generated within CGLS from SPOT- 18 320 VEGETATION (1999-2013) and PROBA-V (2014-2017) satellite imagery. We found that 321 the phenology derived from LAI (or FCOVER) was more closely related to actual ground 322 observation than the NDVI-derived phenology (Table 3). LAI-phenology is based on leaf 323 development rather than on proxies provided by vegetation indices which are not driven solely 324 by the amount of leaves but also by the canopy structure and the leaf biochemical properties 325 (Richardson et al., 2009). LAI is more sensitive than vegetation indices such as NDVI to larger 326 amounts of vegetation (Myneni &Williams, 1994). In addition, NDVI V2.1 is affected by 327 variations in solar zenith and viewing angles and surface reflectance bidirectional effects (Tote 328 et al., 2017). 329 In addition to the vegetation variable, the derived phenology was found to be highly 330 sensitive to the retrieval algorithm and processing chain. We found that the retrieved 331 phenology performed the best using LAI V2 (Table 3) due to the improved continuity (no 332 missing data in V2) and smoothness as compared to V1. In V1 products and in NDVI time- 333 series, the noise primarily due to cloud contamination and atmospheric effects is an important 334 shortcoming in the study of land surface phenology. To overcome this limitation, smoothing 335 methods were applied in SPIRITS approach as a pre-processing step to phenological retrieval. 336 The choice of the smoothing method introduced differences of up to 50% in the performance 337 of the phenology derived from LAI V1 (Table 2). These conclusions agree with previous 338 literature studies highlighting the importance of the temporal reconstruction methods 339 (Kandasamy et al., 2013; Verger et al., 2013). However, the phenology derived from original 340 LAI V2 (smoothing and gap filling was already included in the retrieval algorithm) 341 outperforms the phenology derived from the seasonal trajectories derived after smoothing LAI 342 V1 (Table 3). 343 We tested four state of the art methods to extract phenological metrics: thresholds, logistic 344 function, derivative and moving average. Each method has its own strengths and limitations 345 (de Beurs and Henebry, 2010). The threshold approach based on a percentage of the annual 346 amplitude is simple and robust but it is sensitive to the minimum and maximum values that 347 may be affected by noise in the signal. The logistic function approach has been widely used 348 (e.g. Zhang et al., 2003) but it is limited to the performance of the model fitting (Beck et al., 349 2006) and it may fail when the curvature function is too flat to determine the phenophases (de 350 Beurs and Henebry, 2010). The derivative approach based on the maximum increase and 351 decrease of the vegetation variable is very sensitive to the noise in the signal and the temporal 352 smoothing and composition approach and it cannot represent short growing seasons well, 19 353 especially when the increase and decrease in the annual time series occur rapidly and abruptly 354 (Beck et al., 2006). The moving average approach is based on the assumption that vegetation 355 growth follows a well-defined temporal profile and it may fail in cases of disturbances and 356 abrupt changes. Further, the selection of the time lag is arbitrary. 357 We found that the choice of the extraction method introduced differences >150% and > 358 85% in the performance of the SoS and EoS, respectively, derived from LAI V2 when 359 compared to ground observations (PEP725 and USA-NPN) (Table 4). The percentile method 360 agreed the best with ground measurements. The validation over ground observations indicated 361 that the 30% threshold of the LAI amplitude was optimal for detecting SoS but that a 40% 362 threshold was more suitable for detecting EoS in agreement with Verger et al. (2016). The 363 accuracy for SoS using ground observations produced overall RMSEs of 9 and 11 d for the 364 date of increasing leaf size and leaf unfolding, respectively, for the forests in Europe and the 365 USA. The biases were <2 d (~10 d of standard deviation). We found poorer performances for 366 EoS, than SoS, with higher RMSE of 25 d (28 d) and a bias of 8 d (-6 d) in USA (Europe). 367 The lower performances for the EoS as compared to the SoS is associated to higher 368 uncertainties of both satellite (atmospheric effects, snow and poor illumination conditions) 369 (Delbart et al., 2005) and ground (the timing of leaf colouring is more subjective and difficult 370 to identify than spring phenophases like leaf unfolding) (Estrella & Menzel, 2006) phenology 371 for autumn. Richardson et al. (2009) also reported higher variabilities across the canopy of the 372 timing and rate of foliar development in autumn than spring. 373 The validation of land surface phenology with ground observations presented some 374 difficulties, such as the spatial distribution and the spatial representability of the data. The 375 ground measurements represent the phenology for a limited number of individual plants that 376 are not necessarily the most representative species of the 1-km satellite pixels. Conversely, 377 satellite phenology at 1-km resolution represented an integrated response across landscapes 378 with diverse species and phenological behaviors. The phenology of each species and their 379 characteristics (sizes, ages, homogeneity), though, influenced the satellite signal, depending 380 on its abundance within the pixel sampling area and on the timing of their phenophases 381 (Delbart et al., 2015). Statistics of the comparison between LAI V2 derived EoS using the 382 percentile method and ground measurements improved significantly when the analysis is 383 restricted to Acer rubrum in USA-NPN (RMSE of 10 d, bias of -6 d and significant correlation 384 (p<0.001) of 0.8 (c.f. Table 4)) and to homogeneous sites (forest cover>50% based on GEE 20 385 high resolution imagery) in PEP725 (RMSE of 22 d, bias of 6 d and significant correlation 386 (p<0.05) of 0.4 (c.f. Table 4)). 387 The differences in the definition of the ground phenophases and satellite metrics hamper 388 the comparison. The logistic function and, specially, the moving-average approaches 389 systematically advanced the SoS as compared to ground measurements since these methods 390 determine the SoS as the timing when the vegetation variable starts to increase (Table 4). On 391 the contrary, the derivative approach based on the most rapid increase of the signal introduces 392 positive delays in the SoS. The opposite trend is observed for the EoS (Table 4): the logistic 393 function and moving-average show a delay in the detection of the EoS while the derivatives 394 advances the timing of EoS. On the other hand, ground measurements are subjective and some 395 ambiguity exists in the definition of phenophases. In this sense, we found positive bias of 2 d 396 for the SoS retrieved with the 30th percentile of LAI amplitude when comparing with USA- 397 NPN “leaves” phenophase but negative bias of -2 d when comparing to “increasing leaf size”. 398 Further studies will focus on the comparison of the retrieved land surface phenology with 399 continuous ground observations from PhenoCam (Zhang et al., 2018). This should ultimately 400 lead to propose a standardization in the definition of phenological metrics. 401 This research is expected to contribute for the development of a dedicated algorithm for 402 the operational retrieval of land surface phenology within CGLS. Validation using ground 403 observations was limited to deciduous broadleaf forests. Further studies should extend the 404 analysis to other vegetation types. The methods may need to be adapted to handle multiple 405 and irregularly occurring vegetation growing cycles. Finally, forecasting approaches need to 406 be developed for near-real time land surface phenology retrieval. 407 Acknowledgements 408 This research was supported by an FPU grant (Formación del Profesorado Universitario) 409 from the Spanish Ministry of Education and Professional Training to the first author 410 (FPU2015-497), the Copernicus Global Land Service (CGLOPS-1, 199494-JRC), the Spanish 411 Government grant CGL2016-79835-P, the Catalan Government grant SGR 2017-1005, and 412 the European Research Council Synergy grant ERC-2013-SyG-610028 IMBALANCE-P. 413 414 415 416 21 417 References 418 419 - Atkinson, P.M., Jeganathan, C., Dash, J., Atzberger, C., 2012. Inter-comparison of four 420 models for smoothing satellite sensor time-series data to estimate vegetation phenology. 421 Remote Sens. Environ. 123, 400-417. 422 - Atzberger, C., Eilers, P.H., 2011. Evaluating the effectiveness of smoothing algorithms in 423 the absence of ground reference measurements. Int. J. Remote Sens. 32, 3689-3709. 424 - Baret, F., Weiss, M., Lacaze, R., Camacho, F., Makhmara, H., Pacholcyzk, P., Smets, B., 425 2013. GEOV1: LAI, FAPAR essential climate variables and FCOVER global time series 426 capitalizing over existing products. Part 1: Principles of development and production. Remote 427 Sens. Environ. 137, 299-309. 428 - Beaumont, L., Hartenthaler, T., Keatley, M., Chambers, L., 2015. Shifting time: recent 429 changes to the phenology of Australian species. Clim. Res. 63, 203-214. 430 - Beck, P. S. A., Atzberger, C., Hogda, K. A., Johansen, B., Skidmore, A. K., 2006. Improved 431 monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. 432 Remote Sensing of Environment. 100, 321-334. 433 - Brus, D., Hengeveld, G., Walvoort, D., Goedhart, P., Heidema, A., Nabuurs, G., Gunia, K., 434 2012. Statistical mapping of tree species over Europe. Eur. J. Forest Res. vol. 131, 145-157. 435 - Cai, Z., Jönsson, P., Jin, H., Eklundh, L., 2017. Performance of Smoothing Methods for 436 Reconstructing NDVI Time-Series and Estimating Vegetation Phenology from MODIS Data. 437 Remote Sensing. 9, 1271. 438 - Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple 439 method for reconstructing a high-quality NDVI time-series data set based on the Savitzky– 440 Golay filter. Remote Sens. Environ. 91, 332-344. 441 - de Beurs, K., Henebry, G., 2005. Land surface phenology and temperature variation in the 442 International Geosphere–Biosphere Program high-latitude transects. Global Change Biol. 11, 443 779-790. 22 444 - de Beurs, K. M., Henebry, G. M., 2010. Spatio-temporal statistical methods for modelling 445 land surface phenology. In Hudson, & Keatley (Eds.), Phenological research: Methods for 446 environmental and climate change analysis, (pp. 177–208). 447 - Delbart, N., Kergoat, L., Le Toan, T., Lhermitte, J., Picard, G., 2005. Determination of phe- 448 nological dates in boreal regions using normalized difference water index. Remote Sens. 449 Environ. 97, 26-38. 450 - Delbart, N., Beaubien, E., Kergoat, L., Le Toan, T., 2015. Comparing land surface phenology 451 with leafing and flowering observations from the PlantWatch citizen network. Remote Sens. 452 Environ. 160, 273-280. 453 - Eerens, H., Haesen, D., 2015. User’s manual of Software for the Processing and 454 Interpretation of Remotely Sensed Image Time Series. Version 1.4.0. Available at 455 http://spirits.jrc.ec.europa.eu/files/SpiritsManual_140.pdf 456 - Eerens, H., Haesen, D., Rembold, F., Urbano, F., Tote, C., Bydekerke, L., 2014. Image time 457 series processing for agriculture monitoring. Environ. Model. Soft. 53, 154-162. 458 - Eilers, P. H. C., 2003. A perfect smoother. Analytical Chemistry. 75(14), 3631-3636. 459 - Estrella, N., Menzel, A., 2006. Responses of leaf colouring in four deciduous tree species to 460 climate and weather in Germany. Climate Research. 32, 253-267. 461 - Garrity, S., Bohrer, G., Maurer, K., Mueller, K., Vogel, C., Curtis, P., 2011. A comparison 462 of multiple phenology data sources for estimating seasonal transitions in deciduous forest 463 carbon exchange. Agric. Forest Meteorol. 151, 1741-1752. 464 - Ivits, E., Cherlet, M., Mehl, W., Sommer, S., 2009. Estimating the ecological status and 465 change of riparian zones in Andalusia assessed by multi-temporal AVHHR datasets. Ecol. 466 Indic. 9, 422-431. 467 - Jönsson, P., Eklundh, L., 2002. Seasonality extraction by function fitting to time-series of 468 satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing. 40, 1824-1832. 469 - Kandasamy, S., Baret, F., Verger, A., Neveux, P., Weiss, M., 2013. A comparison of 470 methods for smoothing and gap filling time series of remote sensing observations application 471 to MODIS LAI products. Biogeosciences. 10, 4055-4071. 23 472 - Kandasamy, S., Fernandes, R., 2015. An approach for evaluating the impact of gaps and 473 measurement errors on satellite land surface phenology algorithms: Application to 20year 474 NOAA AVHRR data over Canada. Remote Sensing of Environment. 164, 114-129. 475 - Meier, U., Bleiholder, H., Buhr, L., Feller, C., Hack, H., Heß, M., Lancashire, P.D., Schnock, 476 U., Stauß, R., van den Boom, T., Weber, E., Zwerger, P., 2009. The BBCH system to coding 477 the phenological growth stages of plants history and publications. Journal für Kulturpflanzen. 478 61, 41-52. 479 - Myneni, R. B., Williams, D. L., 1994. On the relationship between FAPAR and NDVI. 480 Remote Sens. Environ. 49, 200-211. 481 - Peñuelas, J., Rutishauser, T., Filella, I., 2009. Phenology feedbacks on climate change. 482 Science. 324 (5929), 887-888. 483 - Reed, B.C., Brown, J.F., Vander, D., 1994. Measuring phenological variability from satellite 484 imagery. J. Veg. Sci. 5, 703-714. 485 - Richardson, A., Braswell, B., Hollinger, D., Jenkins, J., Ollinger, S., 2009. Near-surface 486 remote sensing of spatial and temporal variation in canopy phenology. Ecol. Applic. 19, 1417- 487 1428. 488 - Richardson, A., Keenan, T., Migliavacca, M., Ryu, Y., Sonnentag, O., Toomey, M., 2013. 489 Climate change, phenology, and phenological control of vegetation feedbacks to the climate 490 system. Agric. Forest Meteorol. 169, 156-173. 491 - Schwartz, M., Betancourt, J., Weltzin, J., 2012. From Caprio’s lilacs to the USA National 492 Phenology Network Front. Ecol. Environ. 10, 324-327. 493 - Swets, D. L., Reed, B. C., Rowland, J. R., Marko, S. E., 1999. A weighted least-squares 494 approach to temporal smoothing of NDVI. ASPRS Annual Conference, From Image to 495 Information, Portland, Oregon, May 17-21, Proceedings: Bethesda, Maryland. American 496 Society for Photogrammetry and Remote Sens. 497 - Tateishi, R., Ebata, M., 2004. Analysis of phenological change patterns using 1982-2000 498 Advanced Very High-Resolution Radiometer (AVHRR) data. Remote Sens. 25, 2287-2300. 24 499 - Templ, B., Koch, E., Bolmgren, K., Ungersböck, M., Paul, A., Scheifinger, H., 2018. Pan 500 European Phenological database (PEP725): a single point of access for European data. Int. J. 501 Biometeoroly. 62, 1109-1113. 502 - Toté, C., Swinnen, E., Sterckx, S., Clarijs, D., Quang, C., Maes, R., 2017. Evaluation of the 503 SPOT/VEGETATION collection 3 reprocessed dataset: surface reflectances and NDVI. 504 Remote Sens. Environ. 201, 219-233. 505 - Verger, A., Baret, F., Weiss, M., Kandasamy, S., Vermote, E., 2013. The CACAO method 506 for smoothing, gap filling and characterizing seasonal anomalies in satellite time series. IEEE 507 transactions on Geoscience and Remote Sensing. 51, 1963-1972. 508 - Verger, A., Baret, F., Weiss, M., 2014. Near real time vegetation monitoring at global scale. 509 IEEE J. STARS. 7, 3473-3481. 510 - Verger, A., Bater, M., Weiss, M., Filella., Peñuelas, J., 2015. GEOCLIM A global 511 climatology of LAI, FAPAR, and FCOVER from VEGETATION observations for 1999- 512 2010. Remote Sens. Environ. 166, 126-137. 513 - Verger, A., Filella, I., Baret, F., Peñuelas, J., 2016. Vegetation baseline phenology from 514 kilometric global LAI satellite products. Remote Sens. Environ. 178, 1-14. 515 - Verger, A., Baret, F., Weiss, M., 2017. Algorithm Theoretical Basis Document: LAI, 516 FAPAR, FCOVER Collection 1km, Version 2, Issue I1.40. Available at: https://land. 517 copernicus.eu/global/sites/cgls.vito.be/files/products/GIOGL1_ATBD_FAPAR1kmV2_I1.4 518 0.pdf. 519 - Viovy, N., Arino, O., Belward, A. S., 1992. The best index slope extraction (BISE): A 520 method for reducing noise in NDVI time-series. Int J of Remote Sens. 13, 1585-1590. 521 - White, M. A., de Beurs, K., Didan, K., Inouye, D., Richardson, A., Jensen, O., Lauenroth, 522 W., 2009. Intercomparison, interpretation, and assessment of spring phenology in North 523 America estimated from remote sensing for 1982–2006. Global Change Biol. 15, 2335-2359. 524 - Wu, C., Gonsamo, I., Gough, A., Chen, J., Xu, S., 2014. Modeling growing season 525 phenology in North American forests using seasonal mean vegetation indices from MODIS. 526 Remote Sens. Environ. 147, 79-88. 25 527 - Yu, L., Liu, T., Bu, K., Yan, F., Yang, J., Chang, L., Zhang, S., 2017. Monitoring the long- 528 term vegetation phenology change in northeast China from 1982 to 2015. Scientific Report. 529 7, 14770. 530 - Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F., Gao, F., Huete, A., 531 2003. Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 84, 471-475. 532 - Zhang, X., Friedl, M.A., Schaaf, C.B., 2004. Climate controls on vegetation phenological 533 patterns in northern mid- and high latitudes inferred from MODIS data. Global Change Biol. 534 10, 1133-1145. 535 - Zhang, X., Jayavelu, S., Liu, L., Friedl, M.A., Henebry, G.M., Liu, Y., Schaaf, C.B., 536 Richardson, A.D., Gray, J., 2018. Evaluation of land surface phenology from VIIRS data using 537 time series of PhenoCam imagery. Agr. Forest Meteorol. 256, 137-149. 538 539 26