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

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

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₂ between atmosphere and biosphere [3,4,5,6]. Despite its importance, phenology is not always well described in Earth system models [7,8,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,14,15,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,19,20,21,22,23], which has led to the consolidation of national and continental monitoring networks [24,25,26,27,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₂ 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,38,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 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.

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⁻¹; 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 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.

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₂ and H₂O. Additional vertical CO₂ and H₂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₂ fluxes included (1) coordinate rotation using planar fit method [50]; (2) CO₂ 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₂ was then quality checked [53,54]. The net ecosystem exchange (NEE) flux was corrected by adding storage fluxes (integrated CO₂ fluxes using seven levels of CO₂ profiles when possible, otherwise using 1-point storage) to CO₂ 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²) 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):

$${\mathrm{RGB}}_{\mathrm{DN}}{=\mathrm{R}}_{\mathrm{DN}}+{\mathrm{G}}_{\mathrm{DN}}+{\mathrm{B}}_{\mathrm{DN}}$$

(1)

$$\mathrm{GCC}=\frac{{\mathrm{G}}_{\mathrm{DN}}}{{\mathrm{RGB}}_{\mathrm{DN}}}$$

(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]:

$$\mathrm{CamNDVI}=\frac{{{\mathrm{NIR}}_{\mathrm{DN}}}^{\prime }-{{\mathrm{R}}_{\mathrm{DN}}}^{\prime }}{{{\mathrm{NIR}}_{\mathrm{DN}}}^{\prime }+{{\mathrm{R}}_{\mathrm{DN}}}^{\prime }},$$

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

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

¹ 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²: determination coefficient of linear regression. Regression coefficients are all statistically significant (p < 0.001).

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

$$\mathrm{CamNIRv}=\mathrm{CamNDVI}\times {\mathrm{NIR}}_{\mathrm{DN}}^{\prime },$$

(4)

$$\mathrm{CamRVI}=\frac{{\mathrm{NIR}}_{\mathrm{DN}}^{\prime }}{{\mathrm{R}}_{\mathrm{DN}}^{\prime }},$$

(5)

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.

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:

• We discarded VIs measured with PAR below 600 µmol m⁻² s⁻¹. 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].
• 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 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).
• Daily VIs were gap-filled using the Singular Spectrum Analysis (SSA) method implemented in the “spectral.methods” R package [60].
• 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).

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:

• 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. 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, SOStrs, 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 EOStrs and SOStrs, which is widely used in land surface phenology)

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

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}^{\prime }\right|}\left(\mathrm{day}\right),$$

(6)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}^{\prime }-{\mathrm{y}}_{\mathrm{i}}\right)}^2}}\left(\mathrm{day}\right),$$

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

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 3 and

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

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

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

, 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²) of linear regression between VIs and GPP also range from 0.72 to 0.87, in particular, with slightly higher R² 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. 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) ¹.

	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−2 s−1)	(°C)	(mm)	(µmol m−2 s−1)	(°C)	(mm)	(µmol m−2 s−1)	(°C)	(mm)	(µmol m−2 s−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

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

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.

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

3.3. Comparison of Phneological Transition Dates (PTDs) Derived from VIs and GPP

The comparison of PTDs derived from different VIs and GPP is shown in Figure 6 and Figure 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. 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 ¹.

Season Summary
	GCC		CamNDVI		CamNIRv		CamRVI
Stats	Green-Up	Dry-Down	Green-Up	Dry-Down	Green-Up	Dry-Down	Green-Up	Dry-Down
N	24	30	24	30	24	30	24	30
MAE (day)	9.4	6.3	17.4	4.6	7.0	7.4	20.5	5.0
RMSE (day)	12.8	8.6	20.4	5.5	11.0	8.9	22.7	6.3
p-value	***	***	***	ns	ns	***	***	*

¹ Statistics are computed using the phenological transition dates (PTDs, defined in Table 2) of the Green-up period (UD, SOS_trs, and POS1) and the Dry-down period (EOS_trs, 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.

). 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 (Table 4 and

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

	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

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

).

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

3.4. 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 (Figure 8 and Figure 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).

4. Discussion

4.1. 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 (Figure A2 and Figure 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]).

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 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 (Figure 8, Figure A4 and Figure 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.

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.