Using near-infrared-enabled digital repeat photography to track structural and physiological phenology in mediterranean tree–grass ecosystems

remote sensing Article Using Near-Infrared-Enabled Digital Repeat Photography to Track Structural and Physiological Phenology in Mediterranean Tree–Grass Ecosystems Yunpeng Luo 1,*ID , Tarek S. El-Madany 1ID , Gianluca Filippa 2, Xuanlong Ma 1ID , Bernhard Ahrens 1ID , Arnaud Carrara 3, Rosario Gonzalez-Cascon 4, Edoardo Cremonese 2, Marta Galvagno 2, Tiana W. Hammer 1, Javier Pacheco-Labrador 1,5 ID , M. Pilar Martín5, Gerardo Moreno 6ID , Oscar Perez-Priego 1, Markus Reichstein 1,7, Andrew D. Richardson 8,9,10, Christine Römermann 11,12 ID and Mirco Migliavacca 1 1Department for Biogeochemical Integration, Max-Planck-Institute for Biogeochemistry, 07745 Jena, Germany; [email protected] (T.S.E.-M.); [email protected] (X.M.); [email protected] (B.A.); [email protected] (T.W.H.); [email protected] (J.P.-L.); [email protected] (O.P.-P.); mr[email protected] (M.R.); [email protected] (M.M.) 2Environmental Protection Agency of Aosta Valley, ARPA Valle d’Aosta, 11020 Aosta, Italy; [email protected] (G.F.); [email protected] (E.C.); [email protected] (M.G.) 3Fundación Centro de Estudios Ambientales del Mediterráneo (CEAM), 46980 Paterna, Spain; [email protected] 4 Department of Environment, National Institute for Agriculture and Food Research and Technology (INIA), 28040 Madrid, Spain; [email protected] 5Environmental Remote Sensing and Spectroscopy Laboratory (SpecLab), Spanish National Research Council (CSIC), 28037 Madrid, Spain; mpilar[email protected] 6Institute for Dehesa Research, University of Extremadura, 10600 Plasencia, Spain; [email protected] 7Michael-Stifel-Center Jena for Data-Driven and Simulation Science, 07743 Jena, Germany 8Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, MA 02138, USA; andrew[email protected] 9School of Informatics, Computing and Cyber Systems, Northern Arizona University, Flagstaff, AZ 86011, USA 10 Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ 86011, USA 11 Institute of Ecology and Evolution, Plant Biodiversity Group, Friedrich Schiller University Jena, 07743 Jena, Germany; [email protected] 12 German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, 04103 Leipzig, Germany *Correspondence: [email protected]; Tel.: +49-364-757-6208 Received: 16 July 2018; Accepted: 13 August 2018; Published: 15 August 2018   Abstract: Tree–grass ecosystems are widely distributed. However, their phenology has not yet been fully characterized. The technique of repeated digital photographs for plant phenology monitoring (hereafter referred as PhenoCam) provide opportunities for long-term monitoring of plant phenology, and extracting phenological transition dates (PTDs, e.g., start of the growing season). Here, we aim to evaluate the utility of near-infrared-enabled PhenoCam for monitoring the phenology of structure (i.e., greenness) and physiology (i.e., gross primary productivity—GPP) at four tree–grass Mediterranean sites. We computed four vegetation indexes (VIs) from PhenoCams: (1) green chromatic coordinates (GCC), (2) normalized difference vegetation index (CamNDVI), (3) near-infrared reflectance of vegetation index (CamNIRv), and (4) ratio vegetation index (CamRVI). GPP is derived from eddy covariance flux tower measurement. Then, we extracted PTDs and their uncertainty from different VIs and GPP. The consistency between structural (VIs) and physiological (GPP) phenology was then evaluated. CamNIRv is best at representing the PTDs of GPP during the Green-up period, while CamNDVI is best during the Dry-down period. Moreover, CamNIRv outperforms the other VIs in tracking growing season length of GPP. In summary, the results show it is promising to track structural Remote Sens. 2018,10, 1293; doi:10.3390/rs10081293 www.mdpi.com/journal/remotesensing Remote Sens. 2018,10, 1293 2 of 32 and physiology phenology of seasonally dry Mediterranean ecosystem using near-infrared-enabled PhenoCam. We suggest using multiple VIs to better represent the variation of GPP. Keywords: phenology; tree–grass ecosystem; Dehesa; PhenoCam; near-infrared-enabled digital repeat photography; phenological transition date (PTD); growing season length (GSL) 1. Introduction Phenology is the study of recurring life cycle stages, and their timing and relationship with environmental factors [ 1 , 2 ]. Phenology controls the seasonality of ecosystem functions and plant feedbacks to climate through diverse processes, such as changes in the surface albedo and the exchange of CO 2 between atmosphere and biosphere [ 3 – 6 ]. Despite its importance, phenology is not always well described in Earth system models [ 7 – 9 ], in particular, the environmental factors controlling phenology are still uncertain [ 6 , 10 ]. Therefore, additional efforts to monitor and model plant phenology are needed to improve the representation of phenology in Earth system models [6,11]. Conventional visual monitoring of phenology dates back to 705 CE [ 12 ], and still plays an important role in evaluating the impacts of climate change on ecosystems [ 13 – 16 ]. However, conventional monitoring requires substantial field work, which limits spatial and temporal representativeness [17]. Near-surface remote sensing is becoming a more frequently used tool to monitor vegetation phenology at the ecosystem scale. In recent years, installation of commercial digital cameras for phenology monitoring (i.e., PhenoCam) has proliferated throughout diverse biomes and continents [18–23] , which has led to the consolidation of national and continental monitoring networks [ 24 – 28 ]. The use of PhenoCam consistently reduces manual labor, guarantees time series of high temporal resolution, and creates a permanent data record from which visual interpretation and qualification can be made at any later point in time [ 29 ]. The proximity to the target ecosystem allows the cameras to track phenological transition dates (PTDs), such as leaf emergence, leaf discoloration, senescence, and green up and senescence of vegetation with high temporal resolution [ 23 , 30 ], as well as monitor the different plant types within the camera’s field of view (FOV; [ 24 ]). Nowadays, the increasing number of sites with digital cameras co-located with ecosystem-atmosphere CO 2 flux measurements collected using the eddy covariance (EC) technique are contributing to understanding the relationship between phenology of structure and function of ecosystems [28,31,32]. Green chromatic coordinates (GCC) is the most commonly used vegetation index (VI) extracted from PhenoCam, due to the requirement of only three visible spectral bands for computation, and it is used to represent plant development throughout the season [ 33 ]. PhenoCam, with an additional near-infrared (NIR) band which is more sensitive to vegetation structural change than visible bands, has increased the use of the PhenoCam-based normalized difference vegetation index (CamNDVI) for the same purpose [ 29 , 34 , 35 ]. Both GCC and CamNDVI are considered plausible indexes to bridge satellite and ground-based observations of phenology [ 24 ]. For instance, GCC and CamNDVI have been shown to be effective (though not always consistent) tools for describing greenness variation of individual plant species and ecosystems in a variety of plant functional types [ 29 , 36 ], and for evaluating and linking remote sensing phenology products [ 18 , 37 – 39 ] with ground observations [ 40 ]. Recently, Badgley et al. [ 41 ] introduced a new vegetation index called near-infrared reflectance of vegetation index (NIRv), designed to mitigate the mixed pixel problem (determining the fraction of vegetated land surface and reconstructing the signal attributable to vegetation) to better represent photosynthesis of ecosystems. A strong correlation between satellite-based NIRv and gross primary productivity (GPP) at global scale was observed, which outperforms the correlation between NDVI and GPP [ 40 ]. As such, it would be interesting to know if this new index could provide an advantage in tracking seasonal GPP and phenology compared to the widely used CamNDVI and GCC at Remote Sens. 2018,10, 1293 3 of 32 ecosystem-scale. The computation of NIRv, as well as other VIs, such as the ratio vegetation index (RVI [ 42 ]), is possible with NIR-enabled PhenoCam by following the approach proposed by Petach et al. [ 29 ] and Filippa et al. [ 34 ] for the computation of CamNDVI. However, up to now, we are not aware of studies evaluating the differences between PTDs derived from multiple PhenoCam-based VIs from PhenoCam (GCC, CamNDVI, CamNIRv, and CamRVI). Past studies related to derived PTDs from PhenoCam mainly focus on temperate/boreal forests and grassland (e.g., [ 32 ]), and only a few recent studies have focused on seasonally dry tree–grass ecosystems [ 43 , 44 ]. Considering that tree–grass ecosystems are a widely distributed land cover type, which occupies 16–35% of the Earth’s land surface [ 45 – 47 ], it is necessary to further investigate the methods to extract PTDs for these ecosystems. Moreover, the increasing number of sites with PhenoCam associated with EC flux measurements open interesting perspectives to evaluate: first, the consistency between PTDs derived from PhenoCam-based VIs and PTDs of ecosystem functioning (physiological phenology, i.e., [ 48 ]); second, the direct relationship between PhenoCam-based VIs and GPP. However, to our knowledge, only a few studies pay special attention to the differences between phenology of ecosystem structure and of ecosystem functioning and carbon fluxes [28,32]. In this study, our main objective is to evaluate the potential of PhenoCam to monitor phenology of seasonally dry Mediterranean tree–grass ecosystems. Specifically, the objectives are (1) to characterize structural and physiological phenology of tree–grass ecosystems and their main climatic drivers using PhenoCam and GPP derived from EC measurements; (2) to compare the PTDs and growing season length (GSL) derived from different PhenoCam-based VIs, and to evaluate their performance in tracking the PTDs and GSL derived from GPP. 2. Materials and Methods 2.1. Sites Description, Instrument Set-Up, and Data Sources The sites used in this study are Mediterranean tree–grass ecosystems, composed predominantly of an herbaceous layer and low-density evergreen broadleaf oak trees (Quercus ilex; ~20 tree ha −1 ; Figure 1). Three sites are located approximately 500 m apart from each other in Majadas de Tiétar, Cáceres, Spain (39 ◦ 56’24.68”N, 5 ◦ 46’28.70”W), while one site is located in La Albuera, Spain (38 ◦ 42’6.48”N, 6 ◦ 47’9.24”W). The experimental sites in Majadas de Tiétar belong to a large-scale manipulation experiment, where the three areas of approximately 20 ha were manipulated with addition of nitrogen (FLUXNET ID ES-LM1), nitrogen and phosphorous (FLUXNET ID ES-LM2), and the last was kept as control (FLUXNET ID ES-LMa). The experimental site in the La Albuera (FLUXNET ID ES-Abr) is a natural ecosystem with no manipulation. In this study, we did not focus on the fertilization, but only on the evaluation of the effectiveness of different vegetation indexes to represent the ecosystem functions. The Majadas de Tiétar and La Albuera are characterized by a long-term annual mean air temperature of 16.7 ± 0.2 ◦ C and 18.3 ± 1.5 ◦ C, respectively; while mean annual rainfall is ca. 650 mm and 400 mm, respectively. The rain falls typically from November to May with a very dry summer [49]. In each site, an EC system was installed at 15 m of height to measure the carbon, water and energy fluxes (Section 2.2 for more details). The fluxes data are available from March 2014 in ES-LM1, ES-LM2, and ES-LMa; and from October 2015 for ES-Abr. Two broadband Decagon SRS (spectral reflectance sensor) sensors with a FOV of 36 degrees were installed on a rotating arm in each tower area. Downwelling irradiance and upwelling radiance at 650 nm (red spectral band) and 810 nm (near-infrared spectral band) were measured every 5 min for tree and grasses from 30 October 2015. A NIR-enabled digital camera (Stardot NetCam 5MP), was mounted at the top of the EC tower (facing north) at each site. Images were collected every 30 min (from 10:00 to 14:30 UTC) as JPEG format. The camera settings were defined according to the “PhenoCam” protocol (https://phenocam. sr.unh.edu/webcam/tools/). Sequential red, blue, green (RGB) and RGB + NIR images were collected Remote Sens. 2018,10, 1293 4 of 32 by the Stardot camera according to Petach et al., [ 29 ]. FOVs of cameras in ES-LM1 and ES-LM2 were stable during the study period (from 1 August 2014 to 31 July 2017), whereas the FOV of ES-LMa was not constant and the camera experienced a white balance problem until 3 December 2015. At the ES-Abr site, images were available from 1 January 2016. Hence, RGB and RGB+NIR images were available for the analysis from 1 August 2014 to 31 July 2017 for ES-LM1 and ES-LM2; from 3 December 2015 to 31 July 2017 for the ES-LMa site; and from 1 January 2016 to 31 July 2017 for the ES-Abr site. This guarantees a total of 9 site-years for the following analysis. Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 32 December 2015. At the ES-Abr site, images were available from 1 January 2016. Hence, RGB and RGB+NIR images were available for the analysis from 1 August 2014 to 31 July 2017 for ES-LM1 and ES-LM2; from 3 December 2015 to 31 July 2017 for the ES-LMa site; and from 1 January 2016 to 31 July 2017 for the ES-Abr site. This guarantees a total of 9 site-years for the following analysis. Figure 1. Experimental setup at Majadas de Tiétar and La Albuera in the Spain (left panel) and an example of regions of interest (ROIs) in each experimental site (right panel): Grass, Tree, and Ecosystem ROIs, respectively. At each site, an eddy covariance (EC) system was installed at a height of 15 m to measure the fluxes of the whole ecosystem. A near-infrared-enabled camera was installed at 15 m beside the EC system to take pictures half-hourly between 10:00 and 14:30. Three EC towers are in the Majadas de Tiétar (FLUXNET IDs are ES-LM1, ES-LM2, and ES-LMa, respectively) and an EC tower in the La Albuera (FLUX ID: ES-Abr), respectively (not shown). The map of Majadas de Tiétar was provided courtesy of the Spanish Program of Aerial Orthophotography. 2.2. EC Data Processing and Flux Partitioning to GPP Each EC system consists of a three-dimensional sonic anemometer (R3-50, Gill LTD, Lymington UK) and an infrared gas analyzer (LI-7200, Licor Bioscience, Lincoln, NE, USA) to measure mixing ratios of CO2 and H2O. Additional vertical CO2 and H2O concentration profiles were measured at seven levels between the surface and the measurement height in the EC tower (0.1, 0.5, 1.0, 2.0, 5.0, 9.0, and 15 m above ground with a LI840, Licor Bioscience, Lincoln, NE, USA). Meteorological variables such as air temperature (Ta), wind speed (WS), relative humidity (RH), incoming global radiation (Rg), photosynthetically active radiation (PAR), and precipitation (Prec) were also measured at each site. Raw EC data were collected at 20 Hz, and were then processed using EddyPro 6.2. The main processing procedures for CO2 fluxes included (1) coordinate rotation using planar fit method [50]; (2) CO2 time lag adjustments by covariance maximization in predefined windows; (3) spectral corrections performed for low and high pass-filtering effects according to Moncrieff et al. [51] and Moncrieff et al. [52]. The calculated flux for CO2 was then quality checked [53,54]. The net ecosystem exchange (NEE) flux was corrected by adding storage fluxes (integrated CO2 fluxes using seven levels of CO2 profiles when possible, otherwise using 1-point storage) to CO2 flux. Figure 1. Experimental setup at Majadas de Tiétar and La Albuera in the Spain ( left panel ) and an example of regions of interest (ROIs) in each experimental site ( right panel ): Grass, Tree, and Ecosystem ROIs, respectively. At each site, an eddy covariance (EC) system was installed at a height of 15 m to measure the fluxes of the whole ecosystem. A near-infrared-enabled camera was installed at 15 m beside the EC system to take pictures half-hourly between 10:00 and 14:30. Three EC towers are in the Majadas de Tiétar (FLUXNET IDs are ES-LM1, ES-LM2, and ES-LMa, respectively) and an EC tower in the La Albuera (FLUX ID: ES-Abr), respectively (not shown). The map of Majadas de Tiétar was provided courtesy of the Spanish Program of Aerial Orthophotography. 2.2. EC Data Processing and Flux Partitioning to GPP Each EC system consists of a three-dimensional sonic anemometer (R3-50, Gill LTD, Lymington UK) and an infrared gas analyzer (LI-7200, Licor Bioscience, Lincoln, NE, USA) to measure mixing ratios of CO 2 and H 2 O. Additional vertical CO 2 and H 2 O concentration profiles were measured at seven levels between the surface and the measurement height in the EC tower (0.1, 0.5, 1.0, 2.0, 5.0, 9.0, and 15 m above ground with a LI840, Licor Bioscience, Lincoln, NE, USA). Meteorological variables such as air temperature (Ta), wind speed (WS), relative humidity (RH), incoming global radiation (Rg), photosynthetically active radiation (PAR), and precipitation (Prec) were also measured at each site. Raw EC data were collected at 20 Hz, and were then processed using EddyPro 6.2. The main processing procedures for CO 2 fluxes included (1) coordinate rotation using planar fit method [ 50 ]; Remote Sens. 2018,10, 1293 5 of 32 (2) CO 2 time lag adjustments by covariance maximization in predefined windows; (3) spectral corrections performed for low and high pass-filtering effects according to Moncrieff et al. [ 51 ] and Moncrieff et al. [ 52 ]. The calculated flux for CO 2 was then quality checked [ 53 , 54 ]. The net ecosystem exchange (NEE) flux was corrected by adding storage fluxes (integrated CO 2 fluxes using seven levels of CO2profiles when possible, otherwise using 1-point storage) to CO2flux. The u*-threshold, which was used as a criterion to discriminate lowand well-mixed eddies in the nighttime, was estimated for each year and tower individually (the median of u*-threshold ranges from 0.11 to 0.18 for ES-LM1, ES-LM2, and ES-LMa while 0.20–0.24 for ES-Abr) following Papale et al. [ 55 ]. The time series of NEE were gap filled using the marginal distribution sampling (MDS) method [ 56 ] which is based on lookup tables of temperature, global radiation, and water vapor pressure deficit (VPD) classes for short temporal windows i.e., 14 days. The gap-filled time series of NEE were then partitioned into GPP as described in Reichstein et al. [ 56 ]. In brief, the nighttime flux ( Rg < 10 W/m2) of NEE (i.e., only respiration) is extrapolated from nighttime to daytime through a temperature response function, which is based on short term temperature sensitivities (for details see [ 56 ]). The u*-threshold, gap-filling, and partitioning was performed with the R package REddyproc [57]. 2.3. Calculation of Vegetation Indexes from PhenoCam Digital numbers (DNs) of each individual channel (R DN , G DN , B DN and NIR DN ) were extracted from each photograph, and averaged over the different regions of interest (ROIs) (Figure 1). The overall brightness of each ROI (RGB DN ) and the relative brightness of green channel, known also as green chromatic coordinates—GCC, were computed with Equations (1) and (2): RGBDN=RDN +GDN +BDN (1) GCC =GDN RGBDN (2) CamNDVI was also computed in the different ROIs according to Petach et al. [ 29 ], using the algorithm implemented in the “phenopix” R package [34,36]: CamNDVI =NIRDN0−RDN0 NIRDN0+RDN0, (3) where NIR DN ‘ and R DN ’ are the adjusted exposure values of NIR DN and R DN , respectively. For a detailed calculation and the exposure adjustment formula, please refer to Petach et al. [ 29 ]. As the R DN ’ and NIR DN ’ are not direct measurements of reflectance, the CamNDVI values are not directly comparable to the NDVI from other data sources. Petach et al. [ 29 ] found a linear relationship between CamNDVI and the NDVI derived from the radiometric sensor (ASD FieldSpec 3) using the bands of 750 nm for the NIR and 605 nm for the red. They suggested using the linear regression coefficients to adjust the CamNDVI values for comparability with NDVI from radiometers [29,34]. Therefore, we applied the method suggested by Petach et al. [ 29 ] and Filippa et al. [ 34 ] to rescale CamNDVI using the NDVI derived from the Decagon SRS (VIs was calculated and averaged over a 30 min period to be consistent with VIs from PhenoCam). The coefficients and statistics of the linear regression used to compute the scaling factors are shown in Table 1. In the following only the rescaled CamNDVI values are used and presented. Likewise, the CamNIRv and CamRVI were also calculated with adopting Equations (4) and (5) which refer to Badgley et al. [41] and Chen [42], respectively. CamNIRv =CamNDVI ×NIR0 DN, (4) CamRVI =NIR0 DN R0 DN , (5) Remote Sens. 2018,10, 1293 6 of 32 where NIR DN0 and R DN0 are the adjusted exposure values like Equation (3). A similar approach used for the CamNDVI was used to compute CamNIRv and CamRVI: SRS-based NIRv and RVI were used to adjust the CamNIRv and CamRVI in order to make them comparable with the data derived from other sources (Table 1). The analysis was conducted on various ROIs as depicted in Figure 1: we selected ROIs with only trees, grass, and both (hereafter referred as Tree, Grass, and Eco ROI, respectively). The different sites have different tree/grass proportions in the camera FOVs, as only one direction of ecosystems could be captured from the PhenoCam. However, the fractional tree canopy covers were consistent (~0.20) in the four sites by referring to field surveys and the classification analysis using airborne hyperspectral imagery [ 58 , 59 ]. The analysis of each footprint for each EC tower also indicates GPP is contributed from ~20% tree canopy and ~80% grasses [ 59 ]. In order to reduce the bias introduced by the different ratios of tree/grass in the images and camera FOV, which has to do with logistical constrains during the camera installation, the ecosystem VIs (GCC, CamNDVI, CamNIRv, CamRVI) were computed by using the weighted average of the VIs derived from Grass ROIs of 0.8 and Tree ROIs of 0.2. Table 1. Results of the linear regression between normalized difference vegetation index (NDVI), near-infrared reflectance of vegetation index (NIRv), ratio vegetation index (RVI) retrieved from PhenoCam and spectral reflectance sensor (SRS) at four experimental sites 1. Site VIs N Intercept Slope R2 ES-LM1 NDVI 4000 0.60 1.51 0.69 NIRv 3282 0.06 0.01 0.60 RVI 3389 −2.03 6.61 0.54 ES-LM2 NDVI 4340 0.56 1.12 0.77 NIRv 4223 0.04 0.01 0.55 RVI 4268 −0.14 3.55 0.57 ES-LMa NDVI 3549 0.67 0.89 0.89 NIRv 3522 0.06 0.01 0.73 RVI 3515 −0.20 5.17 0.84 ES-Abr NDVI 3560 0.64 1.18 0.83 NIRv 3361 0.06 0.01 0.82 RVI 3564 −0.78 5.39 0.80 1 ES-LM1, ES-LM2, ES-LMa, and ES-Abr are the FLUXNET-IDs of four experimental sites, respectively. N, number of observations used for each regression; Intercept and slope: the parameters of the linear regression between indexes retrieved from the PhenoCam and SRS; R 2 : determination coefficient of linear regression. Regression coefficients are all statistically significant (p< 0.001). 2.4. Data Filtering and to Compute Daily VIs and GPP After computing the half-hourly VIs (GCC, CamNDVI, CamNIRv, CamRVI) at ecosystem scale, we applied a series of steps to derive robust time series of daily VIs: 1. We discarded VIs measured with PAR below 600 µ mol m −2 s −1 . This procedure was used, on one hand, to filter out the VI values measured during adverse meteorological conditions (i.e., rainy, foggy, or overcast half-hours [ 34 , 48 ]), and on the other hand, Petach et al. [ 29 ] suggests to apply a threshold on PAR to reduce the variability of CamNDVI due to changes in illumination conditions. Here, we selected a more conservative threshold than Petach et al. [29]. 2. Amax.density filter method was developed to filter and retrieve daily VIs. We constructed the probability density function (PDF) of VIs in 3-day moving windows (30 observations), and assigned the value that has highest probability density as the filtered daily value. We did not apply the widely used max method [22,36], which uses the 90th percentile of the VIs value from a 3-day moving window as the filtered daily value. This is because the variability of NIR DN is larger compared to other channels (i.e., R DN , G DN , B DN ) in the PhenoCam, which would result in large variability of VIs (i.e., CamNDVI, CamNIRv and CamRVI) that is especially obvious for Remote Sens. 2018,10, 1293 7 of 32 Mediterranean ecosystems compared to other ecosystems (some comparison using data retrieved from [ 34 ], results are not shown in this study). Hence, we chose to apply the max.density filter to retrieve time series of VIs with less variability, which were not always retrieved by applying max filter in our study. We used an example to demonstrate the better filter performance of max.density compare to max filter methods in our study (Figure A1). 3. Daily VIs were gap-filled using the Singular Spectrum Analysis (SSA) method implemented in the “spectral.methods” R package [60]. 4. Similar to the processing of VIs, the daily GPP was derived from half-hourly data following the step (2) and step (3). 2.5. Phenological Transition Dates (PTDs) Extraction The Mediterranean climate is characterized by rainy late autumn–winters, and warm dry summers. Typically, the studied tree–grass ecosystems are dry and covered with senescent grasses in summer, while they increase in greenness in the late autumn (after the onset of the rainy season) and spring. Considering the characteristics of the phenological cycle described above, we decided to conduct the analysis using the concept of “hydrological years”, which is, here, defined from 1 August to 31 July (Figure 2). Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 32 Figure 2. Conceptual figure of the seasonal dynamic of green chromatic coordinates (GCC), PhenoCam-based normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), ratio vegetation indexes (CamRVI), or daily maximum gross primary productivity (GPP), as well as their phenological transition dates (PTDs) for a “hydrological year”. The black circles represent the original vegetation indexes value (CamGCC, CamNDVI, CamNIRv, CamRVI, or GPP). Results of the smoothing procedure and its uncertainty are shown by red circles and gray area, respectively. The vertical dashed lines represent the PTDs and the corresponding names are shown: UD, SOStrs, POS1, POS2, EOStrs, RD. Two periods are most focused upon in this study: the Green-up period during autumn to winter (green rectangle which including UD, SOStrs, and POS1) and Dry-down period in late spring to summer (light-red rectangle that including POS2, EOStrs and RD). The time interval between UD and RD is defined as the growing season which indicated in the figure. Detailed description of PTDs and the other phenological terms are reported in the Table 2. 2.6. Statistical Analysis All the statistical analyses were conducted with the R 3.4.3 programming language [63]. The differences among PTDs extracted from the different datasets (GCC, CamNDVI, CamNIRv, CamRVI, and GPP) were evaluated using the mean absolute error (MAE) and root mean squared error (RMSE) (Equations (6) and (7)):   n i=1 1' MAE = , i i ? d n ay y -y  (6)     n2 ' i i i=1 1 RMSE= y ? -y , n day  (7) where yi’ and yi were the PTD dates extracted from two different datasets. Wilcoxon signed-rank tests were used to test for statistically significant differences between each paired PTDs from the different datasets given that were not normally distributed, while paired Student’s t tests were used when PTDs were normally distributed. The linear regressions were conducted between time series of VIs and GPP, or between meteorological variables and GPP using ordinary least squares regression (OLS). However, the linear regression between PTDs and GSL extracted from different VIs and GPP was conducted using major axis regression (R package “lmodel2”) to account for errors of similar magnitude in the y and x axis. Figure 2. Conceptual figure of the seasonal dynamic of green chromatic coordinates (GCC), PhenoCam-based normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), ratio vegetation indexes (CamRVI), or daily maximum gross primary productivity (GPP), as well as their phenological transition dates (PTDs) for a “hydrological year”. The black circles represent the original vegetation indexes value (CamGCC, CamNDVI, CamNIRv, CamRVI, or GPP). Results of the smoothing procedure and its uncertainty are shown by red circles and gray area, respectively. The vertical dashed lines represent the PTDs and the corresponding names are shown: UD, SOS trs , POS1, POS2, EOS trs , RD. Two periods are most focused upon in this study: the Green-up period during autumn to winter (green rectangle which including UD, SOS trs , and POS1) and Dry-down period in late spring to summer (light-red rectangle that including POS2, EOS trs and RD). The time interval between UD and RD is defined as the growing season which indicated in the figure. Detailed description of PTDs and the other phenological terms are reported in the Table 2. Remote Sens. 2018,10, 1293 8 of 32 Table 2. Terms used in this study to describe the phenological transition dates (PTDs) and phenological periods. Terms Description Phenological Transition Dates (PTDs) UD Upturn day in the green up period in the autumn SOStrs When 50% of amplitude in the green up period in the autumn is reached POS1 When the first peak of season is reached POS2 When the second peak of season is reached EOStrs When 50% of amplitude in the senescent period in the summer is reached RD Recession day at the end of senescent period in the summer Phenological periods Green-up Greenness/GPP increasing period in the autumn (including UD, SOS trs , and POS1) Dry-down Greenness/GPP decreasing period in the summer (including POS2, EOStrs and RD) GSLRD-UD Growing season length defined in the hydrological year (day length between UD and RD) GSLEOS-SOS Growing season length defined for comparison with GSLRD-UD (day length between EOS trs and SOS trs , which is widely used in land surface phenology) In this study, we developed a PTD extraction method for PhenoCam-based VIs in seasonally dry tree–grass ecosystems. The methodology of PTDs extraction is composed by the following steps: 1. Data were smoothed using the spline method [ 20 , 36 ]; PTDs were extracted using the derivatives of smoothed seasonal cycle [ 61 ] and applying thresholds (i.e., 50%) of amplitude of VIs [ 62 ]. As the start and end of the season are extremely important to characterize the phenology, we defined two sets of PTDs in the start (UD, SOS trs ; Table 2) and end of season (RD, EOS trs ; Table 2) for intercomparison and better characterizing the phenology. These two sets of PTDs are derived based on different perception and methodology. UD and RD are retrieved as the intersection dates between steepest slope and minimum value in the Green-up and Dry-down periods, respectively [ 61 ]. In contrast, SOS trs and EOS trs are retrieved by using the thresholds of 50% amplitude [ 62 ]; i.e., they are defined when 50% of amplitudes are reached in the Green-up and Dry-down periods, respectively. Other extracted PTDs and the phenological periods analyzed in this study were summarized in Figure 2and Table 2. The detailed procedures and corresponding code related to the extraction of PTDs are provided in Appendix B. 2. Uncertainty of extracted PTDs was assessed by extracting PTDs repeatedly (100 times) from an ensemble of time series constructed by summing original data and random noise as described by Filippa et al. [36]. 2.6. Statistical Analysis All the statistical analyses were conducted with the R 3.4.3 programming language [ 63 ]. The differences among PTDs extracted from the different datasets (GCC, CamNDVI, CamNIRv, CamRVI, and GPP) were evaluated using the mean absolute error (MAE) and root mean squared error (RMSE) (Equations (6) and (7)): MAE =1 n n ∑ i=1yi−y0 i(day), (6) RMSE =s1 n n ∑ i=1y0 i−yi2(day), (7) where y i ’ and y i were the PTD dates extracted from two different datasets. Wilcoxon signed-rank tests were used to test for statistically significant differences between each paired PTDs from the different datasets given that were not normally distributed, while paired Student’s ttests were used when PTDs were normally distributed. Remote Sens. 2018,10, 1293 9 of 32 The linear regressions were conducted between time series of VIs and GPP, or between meteorological variables and GPP using ordinary least squares regression (OLS). However, the linear regression between PTDs and GSL extracted from different VIs and GPP was conducted using major axis regression (R package “lmodel2”) to account for errors of similar magnitude in the y and x axis. 3. Results 3.1. Time Series of VIs (GCC, CamNDVI, CamNIRv, CamRVI), GPP, and Their Relationship with Meteorological Conditions The seasonal variation of VIs and GPP, as well as their correlations, are shown in Figure 3and Table 3, respectively. In general, all VIs have distinct seasonal variations and are consistent with the temporal variability of GPP (Figure 3and Table 3). At four sites, all VIs have significant and good correlation with GPP (Pearson’s correlation coefficient: r ≥ 0.85; Table 3). The determination coefficients (R 2 ) of linear regression between VIs and GPP also range from 0.72 to 0.87, in particular, with slightly higher R 2 for CamNDVI and CamRVI. However, VIs have distinct discrepancies on variation range. For instance, the GCC of our sites ranges 0.32–0.42, while the CamRVI could vary between 1 and 7. Distinct interannual variability of VIs and GPP are also observed at the four sites (Figure 3). The time series (GCC, CamNDVI, CamNIRv, CamRVI, and GPP) have an obvious “two-humped” shape in the hydrological year 2014 (Hydro-2014: 1 August 2014–31 July 2015) and hydrological year 2016 (Hydro-2016: 1 August 2016–31 July 2017), while the shape is not clear in the time series for the hydrological year 2015 (Hydro-2015: 1 August 2015–31 July 2016; Figure 3). The absence of the two-humped shape in Hydro-2015 coincides with the warmer mean air temperature (Ta) observed in winter 2015 (9.4 ± 2.4 ◦ C) and lower Ta in spring 2016 (13.5 ± 3.5 ◦ C) compared to Hydro-2014 ( 6.4 ±3.3 ◦C and 16.3 ± 4.3 ◦ C, respectively) and Hydro-2016 (8.0 ± 3.4 ◦ C and 16.4 ±4.6 ◦C , respectively; Table A1 and Figure 4). The onset of growing season in each hydrological year clearly followed the onset of the rainy season, confirming that autumn phenology in these ecosystems is driven by precipitation (Figure 3). Larger precipitation (281.7 mm) in the spring of Hydro-2015 also leads to higher GPP compared to Hydro-2014 (94.9 mm) and Hydro-2016 (94.0 mm; Figure 4and Table A1), both characterized by less than half of the precipitation compared to Hydro-2015. Table 3. Statistics between daily time series of vegetation indexes (PhenoCam-based green chromatic coordinates (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), ratio vegetation index (CamRVI)) and gross primary productivity (GPP) at four experimental sites 1. Site VIs-GPP N r ES-LM1 GCC 1096 0.90 CamNDVI 1096 0.91 CamNIRv 1096 0.93 CamRVI 1096 0.91 ES-LM2 GCC 1096 0.86 CamNDVI 1096 0.87 CamNIRv 1096 0.87 CamRVI 1096 0.87 ES-LMa GCC 607 0.89 CamNDVI 607 0.90 CamNIRv 607 0.85 CamRVI 607 0.90 ES-Abr GCC 635 0.86 CamNDVI 635 0.91 CamNIRv 635 0.86 CamRVI 635 0.91 1 ES-LM1, ES-LM2, ES-LMa, and ES-Abr are the FLUXNET-IDs of four experimental sites, respectively. N, number of daily data; r, Pearson’s correlation coefficients. Remote Sens. 2018,10, 1293 16 of 32 amount in evergreen tree species remain relatively constant throughout the year, and they can utilize their vast root system to access the water deep within the soil [ 64 – 66 ]. By contrast, grasses are highly responsive to rainfall variations in rainy seasons as they tend to use water and nutrients in topsoil with dense shallow roots [ 11 , 67 ]. The yellowing grasses in senescence lose their vigor during the dry and hot summer. Being that the VIs derived from the tree are relatively constant during the year, grasses contribute a large proportion of VIs (Figure A6) and GPP in rainy seasons, but trees contribute more during dry periods [49]. Meteorology plays an important role in governing seasonal variation of VI and GPP ( Figures A2 and A3 ), though the role and importance of water availability and temperature varies across seasons (phenological stages). In autumn, after the dry season, the onset of greenness and GPP is initiated by the onset of the rainy season (Figure 3). In winter, water is not a limiting factor, due to the typical ample precipitation in late autumn and early winter. On the other hand, we observed that winter temperature is an important limiting factor of plant photosynthetic activity [ 68 ] (Figure 4). In spring, with the increase in incoming radiation and day length, temperature, and associated increase in atmospheric evaporative demand (i.e., VPD), we found precipitation strongly correlated with both GCC and GPP (Figure 4); this is consistent with previous findings over Mediterranean ecosystems [ 69 – 71 ]. In summer, which is the dry season, precipitation triggers an abrupt increase of VIs and GPP, as observed after the heavy rain that occurred in the June-2016 and July-2017 (Figure 3). The large rain pulses caused a decrease in temperature and an increase in water availability, further contributing to the regrowth of plants, in particular, at the early stage of the dry season (Figure 3). On the annual scale, growing season lengths (GSLs) derived from ES-Abr tower are significantly shorter than the GSL derived from the other sites in Majadas de Tiétar (i.e., ES-LM1, ES-LM2, and ES-LMa; Figure 8), which might be attributed to a lower water availability at the ES-Abr tower more than the other sites. In fact, ES-Abr is characterized by about 200 mm of rain less than the sites located in Majadas de Tiétar. In this study, we also observed the important influence of both structural (VIs) and physiological phenology (i.e., GPP) exerted by extreme climate events. There was an extremely warm winter followed by wetter than average spring in Hydro-2015 [ 72 ], which led to a significant impact on both VIs and GPP over our study sites (Figure 3and Table A1). The growth and productivity, in particular of the herbaceous layer, was stimulated by high temperature and concomitant water availability in the winter of Hydro-2015, while lower temperature and higher precipitation (~2 times more than Hydro-2014 and Hydro-2016) slow down the growth of plants in the spring (Figure 4and Table A1). This reduced the difference of VIs and GPP between winter and spring, and further caused the observed disappearance of the typical “two-humped” shape of Mediterranean ecosystems (Figure 3). Moreover, the warm spell in winter contributed to an extremely high GPP and greenness (Red dots in Figure 4) that substantially contributed to the annual total GPP [ 72 ]. These results point towards the important, and often neglected, contribution of winter periods to interannual variability of GPP and phenology in the Mediterranean ecosystems. In this study, we did not focus on the effects of the fertilization on PTDs, productivity, and growth; rather, we focused on the development of the framework to characterize phenology in seasonally dry Mediterranean ecosystems using PhenoCams. Therefore, further analysis will focus on better understanding the sensitivity of structural and physiological phenology to nutrient availability, meteorological drivers and rain pulses by means of model–data integration (e.g., [7,20,73,74]). Remote Sens. 2018,10, 1293 17 of 32 4.2. Utilizing Different PhenoCam-Based VIs to Represent Structural Phenology We found distinct differences in the PTDs extracted from GCC and CamNDVI, two widely used indexes from PhenoCam (Figure 5). Overall, the PTDs derived from the GCC anticipate the ones derived from the CamNDVI, which is in line with previous studies [ 28 , 29 , 34 , 75 ]. In the Green-up period (UD, SOS trs , and POS1), changes in CamNDVI are more gradual than the changes in GCC (Figure 5) [ 29 ], which is likely due to the fact that GCC is more sensitive to color changes of leaves under low value of leaf area index (LAI; LAI < 2 [ 28 ]). Wingate et al. [ 28 ] found that the initial increase of GCC during the early growing season is caused by the rapid changes in leaf area and leaf chlorophyll content (Chl). With the increase of foliar biomass associated with shoot elongation and the formation of new leaves [ 28 , 76 ] after the early growing season, CamNDVI, which is more sensitive to the changes of canopy structure, increases continuously while GCC has no substantial changes [ 34 ]. Previous studies also found that GCC becomes saturated when nitrogen content or Chl only reaches half of the maximum value at Green-up periods [ 35 , 75 ]. By contrast, NDVI continues to increase with the maturation of leaves [ 77 ] and the increase of LAI in the canopy [ 35 , 75 ]. This likely explains why there is a big difference in the timing of the first peak (POS1) of greenness between GCC and CamNDVI (Figure 5). Similar to the Green-up period, during the Dry-down period, GCC is more sensitive to color changes in leaves from green to yellow, while CamNDVI is a better proxy for variation of LAI and biomass [ 28 , 29 , 34 ]. During the Dry-down, the total green LAI of grasses start decreasing, but there is a large amount of dry and senescent biomass presented in the top canopy [ 78 ]. As such, we observed an earlier decline of GCC than CamNDVI during the Dry-down period [ 34 ] (Figure 5). The complementary information provided by GCC and CamNDVI point towards the need of using both indices to monitor structural changes of the canopy [ 29 , 34 ], as well as the need for investigating different VIs that can be computed from PhenoCam. To our knowledge, this study is the first attempt to compute the CamNIRv and CamRVI, and to compare them with the more widely used GCC or CamNDVI [ 34 , 41 ]. There is a systematic difference in PTDs extracted from CamNIRv as compared to those derived from other VIs, particularly for PTDs from the Green-up period. The PTDs derived from CamRVI tend to be later than those from CamNDVI and other VIs. This result is in agreement with Viña et al. [ 79 ], which also showed that satellite-based NDVI increases faster than RVI with the increase of LAI in maize and soybean fields. Our results confirmed that PhenoCam-based VIs provide complementary information that can be used to monitor phenology of structure (biomass, greenness). The systematic differences observed between VIs are consistent with results reported in the literature [ 34 , 75 ] and obtained with spectroradiometers [ 29 , 34 ] or satellites [ 80 ]. Future studies are needed to analyze the systematic differences between PhenoCam-based VIs (e.g., the comparison between NIRv and other VIs). For instance, studies that combine physiological measurements and plant traits collected in the field with PhenoCam data [ 75 , 77 ], in parallel with the use of radiative transfer models (e.g., Wingate et al., [ 28 ]) can provide valuable information to better understand the difference between VIs and different aspects of vegetation phenology. 4.3. Combing Different PhenoCam-Based VIs to Represent Physiological Phenology The relationship of VIs derived from PhenoCam imagery with ecosystem-scale carbon fluxes in semi-arid systems is recently receiving more attention. However, to our knowledge, previous studies focused mainly on the relationship between GPP and GCC, and here, we move forward to include other potential VIs that can be derived from PhenoCams. We observed varying performance among PhenoCam-based VIs in tracking physiological phenology as measured by GPP. During Green-up period, the PTDs derived from CamNIRv agreed the best with the PTDs derived from GPP (Figure 6 and Table 4). Apart from CamNIRv, the UD and SOS trs derived from the GCC, CamNDVI, and CamRVI are statistically delayed more than the ones derived from the GPP (Figure 7). Our results about the differences between GPP and GCC are contrasting to previous studies [ 20 , 32 , 75 ], which Remote Sens. 2018,10, 1293 18 of 32 mostly focused on temperate deciduous forests, evergreen needleleaf forests or grasslands, while the relationship between GCC and GPP is comparable to a study focused on a grassland understory of a “open forest savanna” in Australia [ 30 ]. One possible explanation for the discrepancy between previous studies [ 20 , 32 , 75 ] and this study, can be related to the patterns of re-greening of the vegetation in autumn after the onset of the rainy season. At the beginning of the Green-up period, the canopy of the grass is characterized by a relevant amount of dry biomass (~38% senescent grasses as measured between October to November 2014–2015). Therefore, the new emerging grasses contribute to the photosynthetic activities, but not so much to the measureable greenness. For this reason, the GPP increases systematically earlier than GCC. However, there is evidence in the literature that this systematic delay is dependent on the greening patterns and mechanisms of vegetation phenology, therefore, it is vegetation type-dependent [20,32,75]. During the Dry-down period, the PTDs derived from CamNDVI are closer to the ones derived from GPP. The PTDs of EOS trs , and RD derived from GCC and CamNIRv are more advanced than the ones derived from GPP. Dry biomass starts to accumulate in the senescent period at the top of the grass canopy with still a certain amount of living biomass at the bottom. As mentioned above, GCC is sensitive to the changes of color [28,29,34,75], hence, PTDs derived from GCC are anticipated to advance more than the PTDs derived from GPP at Dry-down period. However, more investigations are needed to explain why CamNIRv is also more advanced than GPP at the Dry-Down period. By contrast, PTDs derived from CamNDVI and CamRVI have no statistically significant differences compared to the ones from GPP, which implies the potential to use CamNDVI or CamRVI to represent GPP in the Dry-down period. For the timing of max structural and physiological phenology, we did not observe the systematic differences between PhenoCam-based VIs and GPP for the POS1 (Figure 7). However, the POS1 in our study is not comparable to the timing of the maximum value in other studies, as the POS1 in our study is in the winter time, while others are in the late spring or summer [ 32 , 75 ]. By contrast, the POS2 extracted from GPP is delayed compared to the POS2 extracted from VIs during the spring period (Figure 7), which agrees with previous studies [32,75]. As a key factor controlling net uptake of carbon dioxide [ 81 , 82 ], accurate estimates of GSL has rendered substantial interest, since it has distinct impacts on ecosystem function [ 83 ]. In this study, we also study the GSL as extracted from GPP and PhenoCam-based VIs. GSL derived from CamNIRv is most representative of the GSL derived from GPP (Figures 8,A4 and A5). NIRv is claimed to explain a large fraction of the variance of GPP, and has better representation than NDVI on monthly to annual time scales [ 41 ]. However, from our study, the CamNDVI better tracks the PTDs of GPP more than CamNIRv on the Drying-down period (Figure 7), which makes the GSL derived from the CamNDVI also close to GSL derived from GPP, though with a larger error compared to CamNIRv. Yang et al. [ 35 ] reported high correlation between physiological properties (e.g., leaf nitrogen content, leaf chlorophyll content) and CamNDVI, which implies CamNDVI could potentially track the GPP well. In summary, NIR-enabled PhenoCam-based VIs (e.g., CamNDVI and CamNIRv) can improve the performance of PhenoCam to represent physiological phenology (i.e., variability of PTDs and GSL as derived from GPP). Compared to conventional PhenoCam (only with blue, green, and red channels), NIR-enabled PhenoCam-based VIs take advantage of the fact that green vegetation reflects more NIR than visible light, which makes them more relevant to monitor variation in biomass and seasonal variability in photosynthetic capacity [ 84 ]. More studies investigating other PhenoCam vegetation indices are needed, for instance NIRv and 2-band enhanced vegetation index (EVI2; EVI computed without blue band) could be complementary indexes to be applied to track GPP in future, as both indexes are reported to have a good relationship with GPP [ 41 , 85 , 86 ] and could be computed with two bands (Red and NIR).Our results confirmed that it is promising to utilize the NIR-enabled PhenoCam as a complementary and cost-effective way to characterize GPP, biomass, and phenology. In this study, we mainly focused on seasonal variability (PTDs and GSL) but at shorter time scales; it is still unclear how the different PhenoCam-based VIs presented in this study fully present the variability of GPP. Remote Sens. 2018,10, 1293 19 of 32 For instance, some VIs could not accurately reflect the variation of GPP during short-term changes of weather conditions, like sudden warm spells or droughts, as pointed out in other studies [ 32 , 87 ]. Hence, we strongly suggest using multiple VIs to better characterize the GPP together with additional ancillary measurements, such as meteorology and leaf traits. 5. Conclusions In this study, we assessed the potential to jointly use near-infrared-enabled digital repeat photography and eddy covariance data for monitoring structural and physiological phenology in seasonally dry Mediterranean tree–grass ecosystems. We analyzed 9 site-years using four PhenoCam-based vegetation indices (GCC, CamNDVI, CamNIRv, and CamRVI) and GPP, and we compared the phenological transition dates (PTDs) and growing season length (GSL) derived from the different data streams. We show that, in Mediterranean tree–grass ecosystem, meteorology plays an important role in governing seasonal variation of vegetation indices and GPP, though the importance of water availability and temperature vary across seasons. We show the PTDs derived from VIs differ from each other. For the widely used GCC and CamNDVI, the PTDs extracted from CamNDVI are delayed compared to the ones derived from the GCC, which is likely attributed to GCC, and is more sensitive to color changes, while CamNDVI is more sensitive to LAI and biomass. CamNIRv is best at representing the PTDs of GPP at the Green-up period, while CamNDVI is the best proxy to represent the PTDs of GPP at the Dry-down period. CamNIRv performs best regarding the representation of the GSL of GPP. In summary, we show that it is possible to determine crucial PTDs of structural and physiological phenology through using near-infrared-enabled digital cameras. GPP could be well represented when combining the use of different VIs for this purpose. Author Contributions: Conceptualization, M.M., M.R., C.R., and Y.L.; Methodology, Y.L., G.F., X.M., and M.M.; Formal Analysis, Y.L., and T.S.E.-M.; Supervision, M.M., M.R., C.R., and X.M., Investigation, B.A., R.G.-C., T.W.H., G.M., and M.P.M.; Data Curation, T.S.E.-M., R.G.-C., M.P.M., G.M., and T.W.H.; Writing-Original Draft Preparation, Y.L., Writing-Review & Editing, all coauthors. Acknowledgments: The authors acknowledge the Alexander von Humboldt Foundation for supporting this research with the Max-Planck Prize to Markus Reichstein. Y.L. and M.M. gratefully acknowledge financial support from the China Scholarship Council. We are also thankful for financial support from the Spanish Ministry of Economy and Competitiveness through the FLUXPEC project “Monitoring changes in water and carbon fluxes from remote and proximal sensing in a Mediterranean dehesa ecosystem” (CGL2012-34383). The development of PhenoCam has been supported by the Northeastern States Research Cooperative, NSF’s Macrosystems Biology program (award EF-1065029 and EF-1702697), DOE’s Regional and Global Climate Modeling program (award DE-SC0016011), and the US National Park Service Inventory and Monitoring Program and the USA National Phenology Network (grant number G10AP00129 from the United States Geological Survey). The authors thank Sujan Koirala for relevant technical assistance and comments on graphics. The authors thank two anonymous reviewers and the editor for constructive comments to improve the earlier manuscript. The authors thank Andrew Durso for final proofreading for English. Conflicts of Interest: The authors declare no conflict of interest. Remote Sens. 2018,10, 1293 20 of 32 Appendix A Table A1. Mean air temperature (Ta_mean), sum of precipitation (Prec), and mean gross primary productivity (GPP) of different seasons at Majadas de Tiétar from hydrological year 2014 to 2016 (Hydro-2014 to Hydro-2016) 1. Autumn Winter Spring Summer (September–November) (December–February) (March–May) (June–July) Ta Prec GPP Ta Prec GPP Ta Prec GPP Ta Prec GPP (◦C) (mm) (µmol m−2s−1)(◦C) (mm) (µmol m−2s−1)(◦C) (mm) (µmol m−2s−1)(◦C) (mm) (µmol m−2s−1) Hydro-2014 17.5 333.4 4.7 6.4 105.6 4.3 16.3 94.9 9.4 27.5 67.7 4.2 Hydro-2015 16.8 296.0 3.6 9.4 181.6 5.7 13.5 281.7 11.7 26.4 14.1 7.0 Hydro-2016 17.2 272.7 3.6 8.0 205.9 5.2 16.4 94.0 8.9 27.3 46.4 4.2 1Mean GPP is the average of day time (6:00–18:00) of GPP from Control tower (ES-LMa), Nitrogen tower (ES-LM1), Nitrogen and Phosphorous tower (ES-LM2) at Majadas de Tiétar. Table A2. Comparison between phenological transition dates (PTDs) derived from PhenoCam-based green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), ratio vegetation index (CamRVI) and PTDs derived from GPP in four Mediterranean experimental sites 1. Green-Up Period Dry-Down Period GCC CamNDVI CamNIRv CamRVI GCC CamNDVI CamNIRv CamRVI Stats UD SOStrs UD SOStrs UD SOStrs UD SOStrs EOStrs RD EOStrs RD EOStrs RD EOStrs RD N 8 8 8 8 8 8 8 8 10 10 10 10 10 10 10 10 MAE (day) 8.0 13.1 14 17.5 4.8 2.8 16 20.2 7.0 7.1 4.9 5.5 8.7 8.8 6.9 4.4 RMSE (day) 10.6 10.4 14.5 19 6.3 3.3 16.6 21.6 9.6 10.0 6.0 6.3 10.2 9.9 8.0 5.5 r 0.61 0.64 0.94 0.77 0.91 0.98 0.91 0.75 0.90 0.91 0.94 0.94 0.91 0.95 0.92 0.95 p-value * * *** * ns ns *** ** * * ns ns * * ns ns 1 Statistics are computed using the PTDs at the Green-up period (including PTD UD, SOS trs , and POS1) and the Dry-down period (including PTD EOS trs , RD, and POS2). Here only show the statistics for UD, SOS trs , EOS trs , and RD which are most important in Green-up and Dry-down period. Please refer to each PTD in Table 2. N, number of observations; MAE, mean absolute error; RMSE, root mean squared error. p-values are as follows: *** p≤0.001, ** 0.001 < p≤0.01, * 0.01 < p≤0.05, ns for p> 0.05. Remote Sens. 2018,10, 1293 21 of 32 Remote Sens. 2018, 10, x FOR PEER REVIEW 21 of 32 Figure A1. Comparison between max (Sonnentag et al., 2012) and max.density methods for their performance of filtering the PhenoCam-based vegetation indexes like near-infrared reflectance of vegetation index (NIRv). The black points are original half-hourly NIRv data, while red points stand for the filtered daily NIRv data retrieved by using max or max.density method. Figure A1. Comparison between max (Sonnentag et al., 2012) and max.density methods for their performance of filtering the PhenoCam-based vegetation indexes like near-infrared reflectance of vegetation index (NIRv). The black points are original half-hourly NIRv data, while red points stand for the filtered daily NIRv data retrieved by using max or max.density method. Remote Sens. 2018,10, 1293 22 of 32 Remote Sens. 2018, 10, x FOR PEER REVIEW 22 of 32 Figure A2. Scatter plots between monthly mean air temperature (Ta), monthly sum of precipitation (Prec) and green chromatic coordinate (GCC) during day time of (a) autumn, (b) winter, (c) spring, and (d) summer in hydrological year 2014 to 2016 in four Mediterranean tree–grass ecosystem in Spain. The data in hydrological years of 2014, 2015 and 2016 are colored with green, red and blue color, respectively. The regression line and the formula are shown if the linear regression is statistically significant (p-value: p < 0.05). The gray area represents 95% confidence interval. R2: determination coefficient of linear regression. Figure A2. Scatter plots between monthly mean air temperature (Ta), monthly sum of precipitation (Prec) and green chromatic coordinate (GCC) during day time of ( a ) autumn, ( b ) winter, ( c ) spring, and ( d ) summer in hydrological year 2014 to 2016 in four Mediterranean tree–grass ecosystem in Spain. The data in hydrological years of 2014, 2015 and 2016 are colored with green, red and blue color, respectively. The regression line and the formula are shown if the linear regression is statistically significant (p-value: p< 0.05). The gray area represents 95% confidence interval. R 2 : determination coefficient of linear regression. Remote Sens. 2018,10, 1293 23 of 32 Remote Sens. 2018, 10, x FOR PEER REVIEW 23 of 32 Figure A3. Scatter plots between monthly mean air temperature (Ta), monthly sum of precipitation (Prec) and gross primary productivity (GPP) during day time of (a) autumn, (b) winter, (c) spring, and (d) summer in hydrological year 2014 to 2016. The data in hydrological years of 2014, 2015 and 2016 are colored with green, red and blue color, respectively. The regression line and the formula are shown if the linear regression is statistically significant (p-value: p < 0.05). The gray area represents 95% confidence interval. R2: determination coefficient of linear regression. Figure A3. Scatter plots between monthly mean air temperature (Ta), monthly sum of precipitation (Prec) and gross primary productivity (GPP) during day time of ( a ) autumn, ( b ) winter, ( c ) spring, and ( d ) summer in hydrological year 2014 to 2016. The data in hydrological years of 2014, 2015 and 2016 are colored with green, red and blue color, respectively. The regression line and the formula are shown if the linear regression is statistically significant (p-value: p< 0.05). The gray area represents 95% confidence interval. R2: determination coefficient of linear regression. Remote Sens. 2018,10, 1293 24 of 32 Remote Sens. 2018, 10, x FOR PEER REVIEW 24 of 32 Figure A4. Matrix plots of growing season length (GSLRD-UD) calculated from different data source: gross primary productivity (GPP), PhenoCam-based green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), and ratio vegetation index (CamRVI). Upper right boxes display the scatterplots and 1:1 line. Lower left boxes show the mean absolute error (MAE (days); bold in the center), the Pearson’s correlation coefficient (r; gray colored in the top), and slope of linear regression (slope: gray colored in the bottom) between respective variables. p-values are for correlation test for respective variables and as follows: *** p ≤ 0.001, ** 0.001 < p ≤ 0.01, * 0.01 < p ≤ 0.05. Please see the definition of GSLRD-UD in the Table 2. Figure A4. Matrix plots of growing season length (GSL RD-UD ) calculated from different data source: gross primary productivity (GPP), PhenoCam-based green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), and ratio vegetation index (CamRVI). Upper right boxes display the scatterplots and 1:1 line. Lower left boxes show the mean absolute error (MAE (days); bold in the center), the Pearson’s correlation coefficient (r; gray colored in the top), and slope of linear regression (slope: gray colored in the bottom) between respective variables. p-values are for correlation test for respective variables and as follows: *** p ≤ 0.001, ** 0.001 < p ≤ 0.01, * 0.01 < p ≤ 0.05. Please see the definition of GSL RD-UD in the Table 2. Remote Sens. 2018,10, 1293 25 of 32 Remote Sens. 2018, 10, x FOR PEER REVIEW 25 of 32 Figure A5. Growing season length (GSLEOS-EOS) derived from gross primary productivity (GPP), PhenoCam-based green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), and ratio vegetation index (CamRVI) for four Mediterranean tree–grass ecosystems (FLUXNET IDs are ES-LM1, ES-LM2, ESLMa, and ES-Abr, respectively). The gold squares stand for mean GSL in all site-years, while the black points stand for the GSL derived from ES-Abr. The statistically significant differences were tested between GPP and vegetation indexes—GCC, CamNDVI, CamNIRv and CamRVI. p-values are as follows: ** 0.001 < p ≤ 0.01, * 0.01 < p ≤ 0.05, ns for p > 0.05. Please see the definition of GSLEOS-SOS in the Table 2. Figure A5. Growing season length (GSL EOS-EOS ) derived from gross primary productivity (GPP), PhenoCam-based green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), and ratio vegetation index (CamRVI) for four Mediterranean tree–grass ecosystems (FLUXNET IDs are ES-LM1, ES-LM2, ES-LMa, and ES-Abr, respectively). The gold squares stand for mean GSL in all site-years, while the black points stand for the GSL derived from ES-Abr. The statistically significant differences were tested between GPP and vegetation indexes—GCC, CamNDVI, CamNIRv and CamRVI. p-values are as follows: ** 0.001 < p≤0.01, * 0.01 < p≤0.05, ns for p> 0.05. Please see the definition of GSLEOS-SOS in the Table 2. Remote Sens. 2018,10, 1293 32 of 32 79. Viña, A.; Gitelson, A.A.; Nguy-Robertson, A.L.; Peng, Y. Comparison of different vegetation indices for the remote assessment of green leaf area index of crops. Remote Sens. Environ. 2011,115, 3468–3478. [CrossRef] 80. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O’Connor, B.; Paganini, M. Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island. Remote Sens. Environ. 2018,215, 517–529. [CrossRef] 81. Linderholm, H.W. Growing season changes in the last century. Agric. For. Meteorol. 2006 ,137, 1–14. [CrossRef] 82. Churkina, G.; Schimel, D.; Braswell, B.H.; Xiao, X. Spatial analysis of growing season length control over net ecosystem exchange. Glob. Chang. Boil. 2005,11, 1777–1787. [CrossRef] 83. White, M.; Running, S.; Thornton, P. The impact of growing-season length variability on carbon assimilation and evapotranspiration over 88 years in the eastern US deciduous forest. Int. J. Biometeorol. 1999 ,42, 139–145. [CrossRef] [PubMed] 84. Yengoh, G.T.; Dent, D.; Olsson, L.; Tengberg, A.E.; Tucker, C.J., III. Use of the Normalized Difference Vegetation Index (NDVI) to Assess Land Degradation at Multiple Scales: Current Status, Future Trends, and Practical Considerations; Springer: Berlin, Germany, 2015. 85. Peng, Y.; Gitelson, A.A. Application of chlorophyll-related vegetation indices for remote estimation of maize productivity. Agric. For. Meteorol. 2011,151, 1267–1276. [CrossRef] 86. Jiang, Z.; Huete, A.R.; Didan, K.; Miura, T. Development of a two-band enhanced vegetation index without a blue band. Remote Sens. Environ. 2008,112, 3833–3845. [CrossRef] 87. Yang, J.; Tian, H.; Pan, S.; Chen, G.; Zhang, B.; Dangal, S. Amazon droughts and forest responses: Largely reduced forest photosynthesis but slightly increased canopy greenness during the extreme drought of 2015/2016. Glob. Chang. Boil. 2018. [CrossRef] [PubMed] © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).