Divergent Estimates of Forest Photosynthetic Phenology Using Structural and Physiological Vegetation Indices

The accurate estimation of photosynthetic phenology using vegetation indices (VIs) is important for measuring the interannual variation of atmospheric CO2 concentrations, but the relative performances of structural and physiological VIs remain unclear. We found that structural VIs (normalized difference VI, enhanced VI, and near‐infrared reflectance of vegetation) were suitable for estimating the start of the photosynthetically active season in deciduous broadleaf forests using gross primary production measured by FLUXNET as a benchmark, and a physiological VI (chlorophyll/carotenoid index) was better at identifying the end of the photosynthetically active season for deciduous broadleaf forests and both the start and end of season for evergreen needleleaf forests. The divergent performances were rooted in the combined control of structural and physiological regulations of carbon uptake by plants. Most existing studies of photosynthetic phenology have been based on structural VIs, so we suggest revisiting the dynamics of photosynthetic phenology using physiological VIs, which has significant implications on global plant phenology and carbon uptake studies.


Introduction
Forests account for approximately 4 billion ha (~30%) of the terrestrial surface and offset the abundant anthropogenic CO 2 emissions (Bonan, 2008). The uptake of photosynthetic carbon by forests is strongly seasonal, which substantially influences the interannual variation of atmospheric CO 2 concentrations (Xia et al., 2015). A better understanding of the timing of photosynthetic activity, i.e., photosynthetic phenology, is therefore necessary for improving models of the global terrestrial ecosystem and for more accurately predicting future cycles of climate-carbon feedbacks (Peaucelle et al., 2019;Verger et al., 2015).
The rate of carbon uptake by plants can be represented as the product of absorbed photosynthetically active radiation (APAR) and photosynthetic light-use efficiency (LUE) Penuelas et al., 2011). APAR determines the potential photosynthetic rate and is mainly controlled by canopy structure (Running et al., 2004). The potential photosynthetic rate, however, is often downregulated by environmental stresses, so LUE was introduced to quantify how much of this potential is actually realized (Penuelas et al., 2011). Environmental stresses can induce plant physiological responses, and excess APAR is dissipated by several photoprotective mechanisms, e.g., chlorophyll fluorescence emission and heat dissipation, which are remotely detectable (Demmig-Adams & Adams, 2000;Garbulsky et al., 2014).
Satellite vegetation indices (VIs) have been widely used in recent decades to identify photosynthetic phenology and its response to climate change (Chang et al., 2019;D'Odorico et al., 2015;Gonsamo et al., 2012;Middleton et al., 2016;Wu et al., 2017). The information contents of existing VIs can be categorized into structure and physiology. Structural VIs mainly represent vegetation biomass, with typical examples including the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and near-infrared reflectance of vegetation (NIRv). NDVI and EVI may be the most widely used VIs and are strongly correlated with APAR (Huete et al., 2002). NIRv accounts for the influence of background brightness and quantifies the near-infrared reflectance of terrestrial vegetation (Badgley et al., 2017). NIRv scales well with in situ eddy covariance measurements of CO 2 flux at monthly and longer time scales  because of its sensitivity to structure (Kimm et al., 2020). Physiological VIs such as the photochemical reflectance index (PRI) were designed to characterize plant physiology. PRI represents both the diurnal activity of the xanthophyll cycle and seasonal changes to chlorophyll/carotenoid pigment ratios (Penuelas et al., 1994;Wong & Gamon, 2015b) and has been extensively demonstrated to adequately estimate LUE Penuelas et al., 2011;Zhang et al., 2017).
To the best of our knowledge, most, if not all, studies that have extracted photosynthetic phenology using VIs have relied on structural VIs. These VIs only measure potential photosynthetic rates, overlooking actual rates, so erroneous results are expected. For example, the start of the photosynthetically active season (SOS) for deciduous forests has often been estimated to be later than it was (Jeong et al., 2017;Marien et al., 2019), and both SOS and the end of the photosynthetically active season (EOS) for evergreen forests cannot be accurately estimated using structural VIs (Wu et al., 2017). The main challenge of using physiological VIs, e.g., PRI, to track photosynthetic phenology is the limited availability of the key bands for calculating them. For example, the 531-nm wavelength, which can detect changing pigment levels during the dissipation of heat due to excess APAR (Wong & Gamon, 2015a, 2015b, is specified as an ocean band, and some preprocessing (e.g., atmospheric correction) is not operationally implemented by MODIS products. The calculation of PRI from satellite data sets is therefore time-consuming, hindering its wide application for monitoring photosynthetic phenology.
The recent update of the MODIS reflectance product using the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm has allowed the easy acquisition of surface reflectance for both terrestrial and ocean bands across large areas and at high frequencies of observation (Lyapustin et al., 2012). The advent of this product offers promising prospects for tracking photosynthetic phenology using PRI. MODIS is not equipped with the reference band in the original PRI formula, and alternative bands have been selected to calculate "MODIS PRI" (Goerner et al., 2009;He et al., 2016;Middleton et al., 2016). The "MODIS PRI" calculated from bands 1 and 11 is more closely linked to the seasonal change in chlorophyll/carotenoid pigments than to the daily xanthophyll cycle, and it was recently referred to as the chlorophyll/carotenoid index (CCI) (Gamon et al., 2016).
Structural and physiological VIs can theoretically provide complementary information about photosynthetic activity based on LUE paradigm: the structural VIs represent potential photosynthetic rates, and the physiological VIs characterize the short-term downregulation of the maximum LUE, thus providing a measure of the extent to which the potential is realized under stress (Penuelas et al., 2011). Our understanding of the unique advantages of each VI for estimating photosynthetic phenology and the differences in their performance across different forest types (e.g., deciduous and evergreen) for different phenological indicators (e.g., SOS and EOS), however, is still very limited. To fill this knowledge gap, we assessed and compared the performance of three structural VIs (NDVI, EVI, and NIRv) and one physiological VI (CCI) to extract SOS and EOS for forest photosynthetic activity. All four VIs were calculated using MAIAC reflectance. In situ gross primary productivity (GPP) at the sites of eddy-flux towers in evergreen needleleaf forests (ENFs) and deciduous broadleaf forests (DBFs) was used as a benchmark to validate the VI-derived photosynthetic phenologies.

FLUXNET Data
We selected 33 ENF and 18 DBF sites with eddy-flux towers from the FLUXNET-2015 data set (available at http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/) (Table S1). Detailed information for each site can be found in the supporting information. Sites were selected based on the criteria: (1) they were in the Northern Hemisphere with latitudes >30°(see Figure S1) and (2)  The GPPs in the FLUXNET-2015 data set were calculated using gap-filled data for net ecosystem exchange and the standard flux-partitioning method (Lasslop et al., 2010). We used the GPPs from the nighttime-based approach (Reichstein et al., 2005), where nighttime data were used to parameterize a respiration model that was then applied to the entire data set to estimate ecosystem respiration. GPP was then calculated as the difference between ecosystem respiration and net ecosystem exchange.

MODIS Vegetation Indices
The VIs were calculated using surface reflectance from the MCD19A1 Version 6 product. This product was generated using the MAIAC algorithm, which uses an adaptive time series and spatial analysis to derive atmospheric aerosol concentration and surface reflectance without empirical assumptions (Lyapustin et al., 2012). MCD19A1 provides daily surface reflectance for MODIS bands 1-12 at spatial resolutions of 1 km.
The time series of the MCD19A1 observations were first extracted at the closest pixel to each selected FLUXNET site for comparison with in situ GPP measurements. Pixels contaminated by cloud, snow, or a high aerosol optical depth were then excluded based on the layer of quality assurance of the data set. Data with viewing zenith angles >40°were also excluded to minimize the effects of the bidirectional reflectance distribution function (BRDF) (Middleton et al., 2016;Wang et al., 2020).

Extraction of SOS and EOS
Remotely sensed phenologies using VIs are often biased by the effects of snow (Gonsamo et al., 2012;Jin et al., 2017). Snow generally causes abnormal VI values. All VI values during winter (January, February, and December) were therefore assigned as the 5% percentile of all available high-quality VIs during winters. This "fixed-winter" preprocessing was proposed by (Beck et al., 2006) and can greatly reduce the effects of snow (Miao et al., 2013).
Time series of in situ GPP and the VIs were smoothed using a double logistic function: where O s is the smoothed observation, t is the date (day of year), O b and O e are observations (GPP or VIs) before spring green-up and after senescence, respectively, c and d are the slopes of the first and second inflection points associated with forest growth and recession speed, respectively, p and q are the dates of the two inflection points, and O m is the maximum observation at the peak of growth. The parameters c, d, p, q, and k) were fitted against real observations using least squares.
Two methods were adopted to extract SOS and EOS, and the first is the dynamic-threshold method. In situ smoothed GPP and VIs were normalized using where O n is the normalized observation, O s is the smoothed observation using the double logistic function (Equation 6), and O min and O max are the minimum and maximum smoothed observations, respectively. SOS and EOS were defined as the days of the year when O n crossed a threshold. In this study, we specified the dynamic threshold as 0.2, because it is a commonly used threshold and performs well in the most of cases (Shen et al., 2014).
The method proposed by (Gonsamo et al., 2013) (hereafter referred to as Gonsamo's method) was also employed. Gonsamo's method analytically estimates the SOS and EOS based on the minima and maxima of the first, second and third derivatives, i.e., where p, q, O b , and O m are fitted parameters in Equation 6.
Note that the quality flag of MAIAC is very strict (Lyapustin et al., 2012), which results in low number of high-quality VI observations. We only selected reliable SOS and EOS estimations, such that the number of high-quality observations in a 20-day window around the SOS or EOS are larger than 2 (at least, one before and one after the SOS/EOS). This postprocessing reduced the number of SOS and EOS estimates but assured the reliability of the results.

Results
While with slightly different statistics, SOS and EOS from dynamic-threshold and Gonsamo's methods exhibited very similar patterns. Therefore, for simplicity, we only show here the results from dynamic-threshold method, and the results from the Gonsamo's method can be found in Figures S2 and S3 for DBFs and ENFs, respectively.
SOS and EOS over the DBF sites derived from NDVI, EVI, NIRv, and sCCI were compared with those derived from in situ GPP (see Figure 1). SOS derived from in situ GPP was accurately estimated by all four VIs, albeit with systematically earlier estimates (bias <0). NIRv estimated SOS for DBFs the best, with the highest R 2 (0.70) and the smallest RMSE (9.57 d) and bias (−2.95 d). In contrast, the performance of NDVI was less satisfactory (R 2 = 60, RMSE = 13.14 d and bias = −7.54 d). EVI and sCCI provided intermediate accuracy values with slightly lower R 2 (0.69 and 0.67, respectively) and higher RMSEs (10.98 and 10.37 d, respectively) and biases (−4.95 and −4.65 d, respectively) than for NIRv. The VI-detected EOSs were less correlated with their GPP counterparts (R 2 ≤ 0.27) than were the VI-detected SOS. Estimates of EOS were systematically delayed for all four VIs, but the delay was much shorter for sCCI (3.37 d) than for the other three VIs (≥15.39 d). R 2 (0.27) was higher and RMSE (12.19 d) was smaller for sCCI than the three structural VIs, demonstrating its potential for estimating GPP-derived EOS for the DBFs. Among the three structural VIs, NIRv, and EVI performed comparably, with both much better than NDVI.
Estimates of VI-derived SOS were systematically delayed over the ENF sites ( Figure 2). The physiological VI (i.e., sCCI in this study) resulted in the smallest RMSE (16.61 d) and bias (9.78 d), confirming its overwhelming performance over the three structural VIs. Similarly, the sCCI-based EOS estimates over the ENF sites correlated the best with GPP (R 2 = 0.42) and had the smallest RMSE (17.69 d).

Discussion
The three structural VIs all estimated earlier SOSs over the DBF sites, as also reported by Fu et al. (2014) and Jeong et al. (2017). These earlier estimates were most obvious for NDVI (bias = −7.54 d), due to its high sensitivity to understories. Greening in deciduous forests is often first observed at ground level, e.g., herbs and shrubs, so SOS estimated from NDVI mainly indicates the onset of greening of the understory rather than the trees (Ryu et al., 2014). EVI and NIRv were insensitive to background influences (Badgley et al., 2017;Huete et al., 2002) but had a slightly negative bias (−4.95 and −2.95 d, respectively), because trees require time to increase productivity after foliar unfolding, i.e., the timing of carbon assimilation lags behind structural change in deciduous forests (Kikuzawa, 2003). Contrarily, EOS for the DBF sites was later for the structural VIs compared with GPP, consistent with previous studies (Jeong et al., 2017;Marien et al., 2019). This occurs because photosynthesis in high-latitude northern deciduous forests is mainly stressed by photoperiod, although autumnal leaf fall is affected by temperature variability, causing photosynthesis to end before structural recession (Jeong et al., 2017;Medvigy et al., 2013). Structural VIs are widely used for detecting gradual morphological changes in vegetation, even though ENFs have relatively stable amounts of foliage seasonally. NDVI, EVI, and NIRv are therefore not recommended for detecting the seasonality of GPP in evergreen forests (Gamon et al., 2016;Wang et al., 2020;Wong et al., 2019;Wu et al., 2017) because the derived SOSs and EOSs are highly uncertain (Figure 2).
CCI can detect the physiological dynamics of GPP, unlike the structural VIs. Previous studies have reported that CCI was strongly correlated with the size of the pigment pool (Gamon et al., 2016;Wong et al., 2019), which is an important determiner of LUE and ultimately influences photosynthesis (Croft et al., 2017). The capacity of CCI to track subtle seasonal variations in physiology, specifically the chlorophyll/ carotenoid ratio, was confirmed on both foliar (Wong et al., 2019) and canopy (Gamon et al., 2016;Wang et al., 2020) scales. Our study, for the first time, demonstrated the promising potential of CCI to extract photosynthetic metrics of phenology, especially for those with gradual structural variation, e.g., EOS for DBFs and both SOS and EOS for ENFs.
Forest GPP is simultaneously controlled by its capacities to absorb PAR and to convert APAR into fixed carbon (Penuelas et al., 2011). APAR, representing potential photosynthetic rate, is mainly controlled by forest structure so can be well characterized by structural VIs. Potential photosynthetic rates, however, are often downregulated, which can be tracked by the physiological VIs (Demmig-Adams & Adams, 2000). The structural and physiological controls of the capacity of DBF plants to take up carbon in spring are nearly synchronous, but canopy structural changes are more "visible" from remotely sensed observations (Gamon et al., 2016), so structural VIs are more reliable for estimating SOS for DBFs, especially those insensitive to background influences (e.g., NIRv). The structural recession of DBFs in autumn, however, is gradual, and 10.1029/2020GL089167

Geophysical Research Letters
the photosynthetic rate is mainly controlled by physiology (Gallinat et al., 2015), so EOS for DBFs is detected better by the physiological VIs, e.g., CCI, than by the structural VIs. The canopy structure of evergreens is relatively stable throughout the year, and the dynamics of GPP are mainly determined by physiological constraints (Gamon et al., 2016;Wong & Gamon, 2015b). Physiological VIs may therefore be best for ENFs.
Solar-induced chlorophyll fluorescence (SIF) is also a widely used physiological proxy for tracking GPP dynamics (Jeong et al., 2017;Walther et al., 2016). It is deemed to be directly linked to photosynthetic activity (Porcar-Castell et al., 2014). For example, a recent study demonstrated its high accuracy in capturing EOS (Zhang et al., 2020). Therefore, the comparison between SIF and other VIs is imperative especially for the evergreen species, e.g., tropical forests, for which it is hard to extract their phenology. The temporal and spatial resolutions of satellite-observed SIF are very coarse, which makes the direct comparison between them difficult, while the advent of reconstructed SIF from MODIS data (Zhang et al., 2018) provides a new opportunity to compare the performance between SIF and other VIs.
The structural VIs are also influenced by physiology , and physiological VIs are similarly confounded by canopy structure (Middleton et al., 2016). Pioneering studies have separated the structural and physiological components to obtain more "pure" spectral observations (Hilker et al., 2011;Zeng et al., 2019), which can be potentially used to decouple the structural/physiological influence from physiological/structural VIs. The estimation of photosynthetic phenology is expected to be further improved if such "pure" structural or physiological VIs are adopted.

Conclusions
The capacities of structural (NDVI, EVI, and NIRv) and physiological (CCI) VIs calculated using MAIAC reflectances for estimating forest photosynthetic phenology were assessed and compared over 51 forest sites with eddy-flux towers. Results showed that structural VIs are more suitable for estimating the SOS for DBFs, and the physiology-related CCI is better at identifying the EOS for DBFs and both the SOS and EOS for ENFs. The combined control of structural and physiological regulations on carbon uptake by plants is underlying the divergence.
Our study highlights the unique advantage of CCI over existing structural VIs for estimating photosynthetic phenology, especially for ENFs and EOS for DBFs. Most of the existing trending and attributing studies of photosynthetic phenology have relied on structural VIs. Revisiting global photosynthetic phenology using CCI is therefore necessary, and this will also improve our understanding of the roles of the structural and physiological regulation of global GPP dynamics.