Using Near-Infrared-Enabled Digital Repeat Photography to Track Structural and Physiological Phenology in Mediterranean Tree–Grass Ecosystems

: 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 reﬂectance of vegetation index (CamNIRv), and (4) ratio vegetation index (CamRVI). GPP is derived from eddy covariance ﬂux 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 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.

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][46][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.

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. 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 is 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. The design of the schematic plot of instruments in the experimental sites refers to [17] and the photo of the EC tower inside the graphic was taken by T.S. El-Madany.
The above changes do not affect the scientific results. The manuscript will be updated and the original will remain online on the article webpage, with a reference to this correction. The authors would like to apologize for any inconvenience caused.   [17] and the photo of the EC tower inside the graphic was taken by T.S. El-Madany.

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]; (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 CO 2 profiles when possible, otherwise using 1-point storage) to CO 2 flux.
The u*-threshold, which was used as a criterion to discriminate low-and 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/m 2 ) 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].

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): 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]: 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.
where NIR DN and R DN 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.

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.
A max.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 Remote Sens. 2018, 10, 1293 7 of 32 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 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).

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).   Table 2.

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)): 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. 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 SOS trs 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 EOS trs 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, EOS trs and RD) GSL RD-UD Growing season length defined in the hydrological year (day length between UD and RD)

GSL EOS-SOS
Growing season length defined for comparison with GSL RD-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.
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 2 and 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].

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)): 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 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.

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 3 and Table 3, respectively. In general, all VIs have distinct seasonal variations and are consistent with the temporal variability of GPP ( Figure 3 and 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 (  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 4 and 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 .

Comparison of Phenological Transition Dates (PTDs) Derived from Different VIs
The comparison of PTDs extracted from different VIs is shown in Figure 5. Generally, the PTDs derived from the different VIs show smaller difference between each other at Dry-down periods (POS2, EOStrs, RD) compared to the Green-up periods (UD, SOStrs, POS1). Specifically, PTDs derived from GCC have the smallest differences with PTDs derived from CamNIRv, while they are significantly advanced from the ones from CamNDVI at the Dry-down period ( Figure 5). During Green-up periods, the PTDs derived from GCC are more advanced than the ones from CamNDVI The gray area represents 95% confidence interval. R 2 : determination coefficient of linear regression.

Comparison of Phenological Transition Dates (PTDs) Derived from Different VIs
The comparison of PTDs extracted from different VIs is shown in Figure 5. Generally, the PTDs derived from the different VIs show smaller difference between each other at Dry-down periods (POS2, EOS trs , RD) compared to the Green-up periods (UD, SOS trs , POS1). Specifically, PTDs derived from GCC have the smallest differences with PTDs derived from CamNIRv, while they are significantly advanced from the ones from CamNDVI at the Dry-down period ( Figure 5). During Green-up periods, the PTDs derived from GCC are more advanced than the ones from CamNDVI and CamRVI, while they are also delayed compared to PTDs extracted from CamNIRv. By contrast, the PTDs derived from CamNDVI are all more delayed than the PTDs from CamNIRv, in particular, with more than 10 days difference at Green-up period. However, they have small differences compared to the PTDs extracted from CamRVI (<5 days), but they are significantly advanced at Green-up and delayed at Dry-down periods, respectively. Similar to GCC, the PTDs derived from CamNIRv are also advanced than the ones derived from CamRVI with large difference at Green-up period (>15 days).
Remote Sens. 2018, 10, x FOR PEER REVIEW  12 of 32 and CamRVI, while they are also delayed compared to PTDs extracted from CamNIRv. By contrast, the PTDs derived from CamNDVI are all more delayed than the PTDs from CamNIRv, in particular, with more than 10 days difference at Green-up period. However, they have small differences compared to the PTDs extracted from CamRVI (<5 days), but they are significantly advanced at Green-up and delayed at Dry-down periods, respectively. Similar to GCC, the PTDs derived from CamNIRv are also advanced than the ones derived from CamRVI with large difference at Green-up period (>15 days).  Table 2.

Comparison of Phneological Transition Dates (PTDs) Derived from VIs and GPP
The comparison of PTDs derived from different VIs and GPP is shown in Figures 6 and 7. The PTDs derived from VIs for the Green-up period (UD, SOStrs, and POS1), with the exception of CamNIRv, are systematically delayed from the ones derived from GPP ( Figure 6). By contrast, the

Comparison of Phneological Transition Dates (PTDs) Derived from VIs and GPP
The comparison of PTDs derived from different VIs and GPP is shown in Figures 6 and 7. The PTDs derived from VIs for the Green-up period (UD, SOS trs , and POS1), with the exception of CamNIRv, are systematically delayed from the ones derived from GPP ( Figure 6). By contrast, the PTDs extracted from the VIs have a good agreement with the PTDs extracted from the GPP in the Dry-down period (POS2, EOS trs , RD). The difference between PTDs extracted from VIs and GPP ranges within 10 days in the Dry-down period. By contrast, their difference ranges from 8-25 days for the Green-up period, with the only exception of CamNIRv and GPP, which show a smaller difference (Figure 7 panel UD and EOS trs ).
Specifically, the PTDs derived from CamNIRv at the Green-up period and the ones derived from CamNDVI at the Dry-down period have no significant difference compared to the ones derived from GPP, respectively ( Table 4). The PTDs extracted from CamNIRv and GCC have smaller differences (9.4 and 7.0 days, respectively) with the ones derived from GPP in the Green-up period, while they have larger difference with the PTDs derived from GPP compared to CamNDVI (difference of 4.6 days) and CamRVI (difference of 5.0 days) in the Dry-down period (Tables 4 and A2).
The results show that POS1 derived from VIs have a large difference with the POS1 derived from GPP especially for CamNDVI and CamRVI (>20 days). All the POS2 derived from VIs have no statistically significant difference with the ones extracted from GPP, with the difference between VIs and GPP ranging within 5 days (Figure 7). PTDs extracted from the VIs have a good agreement with the PTDs extracted from the GPP in the Dry-down period (POS2, EOStrs, RD). The difference between PTDs extracted from VIs and GPP ranges within 10 days in the Dry-down period. By contrast, their difference ranges from 8-25 days for the Green-up period, with the only exception of CamNIRv and GPP, which show a smaller difference (Figure 7 panel UD and EOStrs). Specifically, the PTDs derived from CamNIRv at the Green-up period and the ones derived from CamNDVI at the Dry-down period have no significant difference compared to the ones derived from GPP, respectively ( Table 4). The PTDs extracted from CamNIRv and GCC have smaller differences (9.4 and 7.0 days, respectively) with the ones derived from GPP in the Green-up period, while they have larger difference with the PTDs derived from GPP compared to CamNDVI (difference of 4.6 days) and CamRVI (difference of 5.0 days) in the Dry-down period (Tables 4 and A2).
The results show that POS1 derived from VIs have a large difference with the POS1 derived from GPP especially for CamNDVI and CamRVI (>20 days). All the POS2 derived from VIs have no statistically significant difference with the ones extracted from GPP, with the difference between VIs and GPP ranging within 5 days (Figure 7).    Table 2) of the Green-up period (UD, SOStrs, and POS1) and the Dry-down period (EOStrs, RD, and POS2). N, number of observations; MAE, mean absolute error; RMSE, root mean squared error. p-values are as follows: *** p ≤ 0.001, * 0.01 < p ≤ 0.05, ns for p > 0.05.  Table 2.

Comparison of Growing Season Length (GSL)-Derived VIs and GPP
The growing season length (GSL: defined as the difference between PTDs of RD and UD) derived from GPP shows no statistically significant difference with GSL calculated from CamNIRv ( Figure 8). By contrast, GSLs derived from GCC, CamNDVI, and CamRVI is statistically significantly different from GSL derived from GPP with mean absolute error (MAE) of 12.2, 10, and 15.5 days, respectively ( Figure A4). We also compare the GSLs derived from GPP and VIs with the GSL defined as the difference between EOStrs and SOStrs, which is widely used in the remote sensing field. The results ( Figure A5) also agree with the above results.

Figure 7. Barplots of the differences between phenological transition dates (PTDs) extracted from different PhenoCam-based vegetation indexes (green chromatic coordinate (GCC), normalized difference vegetation index (CamNDVI), near-infrared reflectance of vegetation index (CamNIRv), ratio vegetation index (CamRVI)
) and gross primary productivity (GPP). The error bars represent the standard error of the PTDs derived from different experimental sites, and hydrological years. The statistically significant differences were tested using Wilcoxon signed-rank tests (when PTDs were not normally distributed) and paired Student's t tests (when PTDs were normally distributed). p-values are as follows: *** p ≤ 0.001, ** 0.001 < p ≤ 0.01, * 0.01 < p ≤ 0.05, ns for p > 0.05. Please see the definition of each PTD in the Table 2. Table 4. 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 .

Comparison of Growing Season Length (GSL)-Derived VIs and GPP
The growing season length (GSL: defined as the difference between PTDs of RD and UD) derived from GPP shows no statistically significant difference with GSL calculated from CamNIRv ( Figure 8). By contrast, GSLs derived from GCC, CamNDVI, and CamRVI is statistically significantly different from GSL derived from GPP with mean absolute error (MAE) of 12.2, 10, and 15.5 days, respectively ( Figure A4). We also compare the GSLs derived from GPP and VIs with the GSL defined as the difference between EOS trs and SOS trs , which is widely used in the remote sensing field. The results ( Figure A5) also agree with the above results. In general, the GSL calculated from GPP is larger than the ones derived from the VIs (Figures 8  and A4). Most of GSLs range from 240 to 270 days, with the exception of GSL derived from ES-Abr sites, which is the driest site among the four studied ( Figure A4). GSLs derived from different VIs have good correlation between each other (r >0.8; Figure A4). Specifically, GSL derived from GCC has the least difference (MAE: 4.8 days) between the ones derived from CamNDVI. GSL derived from CamNDVI and CamRVI has a difference of 5.8 days. Contrastively, GSL derived from GCC, CamNDVI and CamRVI has a relatively high difference (>8 days) between the GSL derived from the CamNIRv ( Figure A4).  Table 2.

Characterizing Variatoin and Drivers of Structural and Physiological Phenology
Here, we discuss and characterize the phenology of a seasonally dry Mediterranean tree-grass ecosystems using high temporal resolution of VIs (i.e., daily) derived from PhenoCam and GPP from EC towers. We found large seasonal variations in both PhenoCam-based VIs (structure) and GPP (physiology). In general, seasonal variations in PhenoCam VIs were in phase with those of GPP ( Figure 3 and Table 3). We argue that seasonal variations in VIs and GPP are driven by the herbaceous layer, which dominates the ecosystem dynamics in our study sites [59]. Our sites are characterized by relatively sparse evergreen broadleaf trees (~20%) and a larger fraction of annual grasses [59]. Foliage 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][65][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]. In general, the GSL calculated from GPP is larger than the ones derived from the VIs (Figures 8 and A4). Most of GSLs range from 240 to 270 days, with the exception of GSL derived from ES-Abr sites, which is the driest site among the four studied ( Figure A4). GSLs derived from different VIs have good correlation between each other (r >0.8; Figure A4). Specifically, GSL derived from GCC has the least difference (MAE: 4.8 days) between the ones derived from CamNDVI. GSL derived from CamNDVI and CamRVI has a difference of 5.8 days. Contrastively, GSL derived from GCC, CamNDVI and CamRVI has a relatively high difference (>8 days) between the GSL derived from the CamNIRv ( Figure A4).

Characterizing Variatoin and Drivers of Structural and Physiological Phenology
Here, we discuss and characterize the phenology of a seasonally dry Mediterranean tree-grass ecosystems using high temporal resolution of VIs (i.e., daily) derived from PhenoCam and GPP from EC towers. We found large seasonal variations in both PhenoCam-based VIs (structure) and GPP (physiology). In general, seasonal variations in PhenoCam VIs were in phase with those of GPP ( Figure 3 and Table 3). We argue that seasonal variations in VIs and GPP are driven by the herbaceous layer, which dominates the ecosystem dynamics in our study sites [59]. Our sites are characterized by relatively sparse evergreen broadleaf trees (~20%) and a larger fraction of annual grasses [59]. Foliage 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][65][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][70][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 3 and 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 4 and 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]).

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.

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 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.
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.

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.           The procedures to extract PTDs from time series of vegetation indexes (VIs) or gross primary productivity (GPP) are as follows (also referring to Figure B1:


Try to find interminD. The linking point of two peaks (interminD) is the minimum between the two peaks between the Julian day of the year (Doy) 150-250.  Once we found the interminD, the time series is split into two parts and for each part the main PTDs are computed. The date of POS is calculated by determining the date at which the maximum value of the time series (using f(t) to refer to the time series hereafter) is reached. Baseline and maxline are the minimum and maximum value in each part of f(t), respectively.  The maximum and minimum of the first derivative of the f(t) (f'(t)) represent the maximum slopes of the upward and downward period (dashed lines). The intersections between these lines and the baseline are defined as upward day (UD) and recession day (RD). UD stands for the value when the f(t) begins to increase during the Green-up period. RD stands for the value when the f(t) stops decreasing during the Dry-down period. The intersections between these lines and maxline are the saturating day (SD) and downward day (DD). The SD indicates when the plants begin to reach full greenness or maximum photosynthesis, while DD stands for the date when plants begin to senesce [61].  SOStrs and EOStrs are retrieved by computing the date when the value reaches 50% of the maximum in the upward and downward period, respectively [62].

Appendix B. The Procedures and Code for Extracting Phenological Transition Dates (PTDs)
The procedures to extract PTDs from time series of vegetation indexes (VIs) or gross primary productivity (GPP) are as follows (also referring to Figure A7:

•
Try to find interminD. The linking point of two peaks (interminD) is the minimum between the two peaks between the Julian day of the year (Doy) 150-250. • Once we found the interminD, the time series is split into two parts and for each part the main PTDs are computed. The date of POS is calculated by determining the date at which the maximum value of the time series (using f(t) to refer to the time series hereafter) is reached. Baseline and maxline are the minimum and maximum value in each part of f(t), respectively.

•
The maximum and minimum of the first derivative of the f(t) (f'(t)) represent the maximum slopes of the upward and downward period (dashed lines). The intersections between these lines and the baseline are defined as upward day (UD) and recession day (RD). UD stands for the value when the f(t) begins to increase during the Green-up period. RD stands for the value when the f(t) stops decreasing during the Dry-down period. The intersections between these lines and maxline are the saturating day (SD) and downward day (DD). The SD indicates when the plants begin to reach full greenness or maximum photosynthesis, while DD stands for the date when plants begin to senesce [61].
• SOS trs and EOS trs are retrieved by computing the date when the value reaches 50% of the maximum in the upward and downward period, respectively [62].
The Code for extracting PTDs is shared in the attached file (PTD_extraction.R). You can also use the example in the attached example folder for a test to extract PTDs.
Remote Sens. 2018, 10, x FOR PEER REVIEW 27 of 32 The Code for extracting PTDs is shared in the attached file (PTD_extraction.R). You can also use the example in the attached example folder for a test to extract PTDs.   Table 2, respectively.