An NDVI-Based Vegetation Phenology Is Improved to be More Consistent with Photosynthesis Dynamics through Applying a Light Use Efficiency Model over Boreal High-Latitude Forests

Remote sensing of high-latitude forests phenology is essential for understanding the global carbon cycle and the response of vegetation to climate change. The normalized difference vegetation index (NDVI) has long been used to study boreal evergreen needleleaf forests (ENF) and deciduous broadleaf forests. However, the NDVI-based growing season is generally reported to be longer than that based on gross primary production (GPP), which can be attributed to the difference between greenness and photosynthesis. Instead of introducing environmental factors such as land surface or air temperature like previous studies, this study attempts to make VI-based phenology more consistent with photosynthesis dynamics through applying a light use efficiency model. NDVI (MOD13C2) was used as a proxy for both fractional of absorbed photosynthetically active radiation (APAR) and light use efficiency at seasonal time scale. Results show that VI-based phenology is improved towards tracking seasonal GPP changes more precisely after applying the light use efficiency model compared to raw NDVI or APAR, especially over ENF.


Introduction
Phenology is the study of the timing of vegetation growing/senescence events based on vegetation seasonal dynamics [1]. Satellite observations of land surface provides spatially-continuous time-resolved reflectance data of terrestrial vegetation. Therefore, extracting phenological timing from time-series remote sensing data is an effective tool for monitoring large-scale vegetation dynamics, especially where in-situ observations are lacking (e.g., [2]). Generally, remote sensing of vegetation phenology is based on the seasonal cycle of reflectance-based vegetation indices (VIs). VIs have been regarded as proxies of the green biomass of vegetation canopy, carrying information of leaf greenness and canopy structure [3]. Because VIs often show clear seasonal cycles over mid-high latitude forests, due to seasonal changes of temperature and precipitation [4,5], numerous researches have utilized VI time-series curves to study the phenology of such vegetation.
Among all the reflectance-based VIs, the normalized difference vegetation index (NDVI) has been most widely used to understand vegetation phenology at both regional and global scales. NDVI is a proxy of green biomass, which is related to canopy photosynthesis. Many studies used NDVI to track seasonal dynamics of carbon assimilation, comparing NDVI with gross primary production (GPP) which represents the photosynthetic carbon uptake [6,7] over boreal forests. However, NDVI-based durations of growing season and photosynthesis season have been found to be significantly longer than those derived from both tower-based and simulated GPP [8,9]. As a matter of fact, NDVI increases earlier in the spring and decreases later in the fall than GPP. As a result, start of season (SOS) date determined by NDVI is earlier than that determined by GPP while end of season (EOS) date appears to be later. Thus, NDVI is considered not sensitive for tracking photosynthesis dynamics of vegetation alone.
However, a direct comparison between NDVI-based phenology and GPP-based phenology may be not reasonable. GPP can be described by the light use efficiency (LUE) model as follows: where PAR stands for photosynthetically active radiation (i.e., shortwave solar radiation between 400 nm and 700 nm) and fPAR stands for fraction of absorbed PAR (APAR). LUE model is a conceptual model originated from [10] based on physiological process of light absorption (APAR) and conversion (LUE) [11]. According to the LUE model, dynamics of GPP are jointly controlled by changes of PAR, fPAR and LUE. At seasonal scale, PAR is related to solar zenith angle and fPAR is determined by canopy structure and leaf greenness, which are directly related to NDVI, while LUE dynamics are more complex because it is related not only to vegetation growing status but to environmental (e.g., radiation, temperature and soil water) stress as well. Thus, in order to use NDVI-based phenology to track photosynthesis dynamics, one should also take seasonal cycles of PAR and LUE into consideration. Because VIs have been used as a robust proxy of fPAR [12] while PAR can be obtained from external data sources, many studies have used VI × PAR as a proxy for APAR [8,9]. However, Moderate Resolution Imaging Spectroradiometer-(MODIS-, onboard Aqua and Terra satellite instruments) based APAR is found to overestimate both the greenness season duration and the photosynthesis duration [8], and APAR is temporally biased compared to GPP [13], indicating that VI-based APAR is still not an appropriate proxy for carbon assimilation. Thus, LUE might be the main reason for inconsistencies between APAR-based phenology and GPP-based phenology.
LUE can be estimated through narrow-band VIs such as photochemical reflectance index (PRI) [14,15] or by introducing maximum LUE, temperature and precipitation stress factors [16]. However, seasonal changes of LUE that are similar to seasonal VI patterns have been observed by studies over northern hemisphere plants (e.g., [17]). Although directly estimating LUE from VI lacks physically-based theoretical basis, some empirical relation may exist between seasonal VI and LUE. Stagakis, et al. [18] reported a significant linear correlation (determination coefficient, r 2 = 0.78) between seasonal NDVI and LUE. Wu et al. [19] found a significant linear relationship between VIs (NDVI and enhanced vegetation index, EVI) and LUE over maize and then used VIs as proxies for LUE at seasonal scale. Based on this scope, a VI-based LUE model is developed: where VI could be either NDVI or EVI and m is an empirical adjust parameter which is determined from the calibration dataset. Favorable results have been reported for this model by several studies (e.g., [19,20]). Most of the applications of this model on GPP estimation are site-based, where in-situ flux-based GPP data that are used to calculate parameter m are available. However, for phenology applications, it is not necessary to calculate m from calibration datasets because in these situations we focus on the seasonal fluctuations of the indicators we use rather than their absolute values. If VI is able to be a proxy of LUE at seasonal scale, as is the case reported by previous studies, the VI-based LUE model should be able to improve the VI-based phenology towards tracking photosynthesis dynamics more precisely. Based on this assumption, this study aims to evaluate the performance of VI-based LUE model in tracking the seasonal cycle of GPP and to use this model to improve VI-based vegetation phenology.

Region and Vegetation Types of Interest
In this study, we focused on the high-latitude forests (between 50 • N and 60 • N) over the Northern Hemisphere. This area mainly includes evergreen forests and deciduous forests. These forests play an important role in the global carbon cycle [21], making a substantial contribution to global carbon sink. Besides, these forests are sensitive to global warming, showing an extending growing season under the influence of increasing winter and the spring temperatures [22]. The forests are divided into three sets according to their region and vegetation type, i.e., North American evergreen needleleaf forests (ENF), Europe-Asia (Eurasia) ENF and Eurasia deciduous broadleaf forests (DBF). The landcover information was obtained from MODIS land cover product (MCD12C1, 2011).

NDVI, GPP, and Radiation Datasets
MODIS vegetation index product (MOD13C2) was used in this study because compared to other datasets (e.g., GIMMS NDVI3g), MODIS NDVI is based on narrow-band calculations specifically designed for vegetation monitoring. In addition, MOD13C2 is quality controlled global-gridded (0.05 • × 0.05 • spatial resolution) monthly average product, which benefits from improved radiometric sensitivity and atmospheric/geometric (bi-directional radiation distribution function, BRDF) corrections [23].
For GPP, we used monthly global products provided by the Max Planck Institute's Biogeochemical integration group (MPI-BGC). The MPI-BGC GPP products were generated through a machine learning algorithm upscaling flux-tower-based GPP into global 0.5 • × 0.5 • products [24]. Here we did not use MODIS GPP products (MOD17A2) because MOD17A2 is also based on LUE model and takes MODIS NDVI as a major input. Thus results might be biased when comparing MODIS GPP and VI-based LUE model.
The incoming shortwave radiation (SWR) was used in this study as a proxy for PAR. Monthly global SWR data (gridded 1 • × 1 • ) provided by the National Aeronautics and Space Administration Clouds and the Earth's Radiant Energy System (NASA-CERES, a series of satellite instruments) are available at https://eosweb.larc.nasa.gov/order-data.

Data Processing and Growing Season Determination
Monthly NDVI, GPP and radiation data covering a three-year period from 2009 to 2011 were used in this study. The temporal range was determined according to time (the end of 2011) of the latest publicly-available MPI-BGC GPP product. VI and radiation products were resampled to the spatial resolution of GPP product. MOD13C2 products were resampled by mean of regional averaging and CERES SWR products were resample by mean of bilinear interpolation. Then single-pixel-based temporal analyses were conducted at the spatial resolution of 0.
To determine the growing season, we used a method proposed by Jeong et al. [22]. This method calculates the dates for SOS and EOS using a threshold of annual maximum of the selected indicators which are NDVI, NDVI × PAR, NDVI 2 × PAR and GPP for this study. Seasonal values were firstly normalized to the range [0, 1] using annual maximum and minimum, then a 6-degree polynomial equation was used to fit the seasonal cycle curve and finally the spring and fall dates were calculated according to the selected threshold. In this study we used a threshold of 0.5.

Results
Four indicators were used in this study to represent the seasonal vegetation dynamics, including NDVI, NDVI × PAR, NDVI 2 × PAR and GPP. Normalized regional-mean seasonal cycles from 2009 to 2011 of the selected indicators are shown in Figure 1 by different regions and vegetation types. According to Figure 1, seasonal patterns of the four indicators stayed similar across years, showing a clear seasonal cycle with single-peak around July. North America ENF and Eurasia ENF showed similar patterns while results for Eurasia DBF show some differences. For North America and Eurasia ENF, the NDVI 2 × PAR curves nearly perfectly matched the GPP curves, especially during the growing and senescence season. The small discrepancies between the two kinds of curves over peak growing season (June-August) might result from the saturation effect of NDVI. Both NDVI and NDVI × PAR curves were found to overestimate the length of growing season, i.e., inclined earlier while declined later than NDVI 2 × PAR and GPP. More discrepancies were found in the fall than in the spring among the seasonal patterns. Much steeper spring increases than fall decreases were found for NDVI over years while the other three indicators showed relatively symmetrical growing and senescence patterns (Figure 1a-f) indicating that NDVI shows more significant lag to photosynthesis (GPP) changes in the fall. Although NDVI × PAR (as a proxy for APAR) curves were much closer to NDVI 2 × PAR and GPP curves in the fall, which was consistent with previous study [9], APAR failed to alleviate but tended to aggravate the early increase in the spring (Figure 1a,c,d). Similar results were reported by Walther et al. [13] when using EVI instead of NDVI. In contrast, NDVI 2 × PAR curves correspond well to GPP curves in both spring and fall, indicating that besides APAR, LUE also plays an important role in the spring recovery of photosynthesis and the spring LUE changes might be tracked with NDVI. The similar patterns in the fall for the three indicators except for NDVI also indicate that radiation reduction is the main factor leading photosynthesis decrease during that season. Thus, the late fall decrease of NDVI compared to GPP, which is often reported (e.g., [2]), is expected to be alleviated when introducing radiation constraint. In addition, the regional standard deviations of NDVI are largely reduced after multiplying with radiation, indicating that more robust According to Figure 1, seasonal patterns of the four indicators stayed similar across years, showing a clear seasonal cycle with single-peak around July. North America ENF and Eurasia ENF showed similar patterns while results for Eurasia DBF show some differences. For North America and Eurasia ENF, the NDVI 2 × PAR curves nearly perfectly matched the GPP curves, especially during the growing and senescence season. The small discrepancies between the two kinds of curves over peak growing season (June-August) might result from the saturation effect of NDVI. Both NDVI and NDVI × PAR curves were found to overestimate the length of growing season, i.e., inclined earlier while declined later than NDVI 2 × PAR and GPP. More discrepancies were found in the fall than in the spring among the seasonal patterns. Much steeper spring increases than fall decreases were found for NDVI over years while the other three indicators showed relatively symmetrical growing and senescence patterns (Figure 1a-f) indicating that NDVI shows more significant lag to photosynthesis (GPP) changes in the fall. Although NDVI × PAR (as a proxy for APAR) curves were much closer to NDVI 2 × PAR and GPP curves in the fall, which was consistent with previous study [9], APAR failed to alleviate but tended to aggravate the early increase in the spring (Figure 1a,c,d). Similar results were reported by Walther et al. [13] when using EVI instead of NDVI. In contrast, NDVI 2 × PAR curves correspond well to GPP curves in both spring and fall, indicating that besides APAR, LUE also plays an important role in the spring recovery of photosynthesis and the spring LUE changes might be tracked with NDVI. The similar patterns in the fall for the three indicators except for NDVI also indicate that radiation reduction is the main factor leading photosynthesis decrease during that season. Thus, the late fall decrease of NDVI compared to GPP, which is often reported (e.g., [2]), is expected to be alleviated when introducing radiation constraint. In addition, the regional standard deviations of NDVI are largely reduced after multiplying with radiation, indicating that more robust results are expected when using APAR or VI-based LUE model to monitor forests' seasonal dynamics rather than using NDVI itself.
Compared to ENF, DBF shows narrower seasonal cycles (Figure 1g-i). This is due to the fact that in the spring DBF rely on newly grown leaves for photosynthesis while new leaves only give limited contribution to ENF's total photosynthesis, while in the fall DBF experiences more sharp decline of canopy greenness than ENF. Both NDVI and APAR show broader seasonal cycles than GPP. It is shown that for DBF, the NDVI 2 × PAR curves still slightly overestimate the duration of growing season, especially in the spring. This is possibly because photosynthesis for DBF starts only when greenness and temperature increase to a higher level than ENF, leading more lag between spring VIs and GPP. It might be interpreted that it is more difficult to use VIs to track seasonal LUE changes in DBF. One possible explanation lies in that for DBF, seasonal variations of NDVI are jointly affected by changes of canopy structure and leaf pigment while for ENF whose canopy structure shows less seasonal dependence, NDVI variations are mainly explained by changes of canopy pigment concentration, which is directly related to LUE. In this sense, NDVI is more reasonable a proxy for LUE in ENF, where the VI-based LUE model performs better in the present study, than in DBF. Although monitoring phenology over ENF is generally more difficult than over DBF using remote sensing (e.g., [25]), the results in the present study suggest that the VI-based LUE model can track GPP dynamics more precisely in ENF.
SOS and EOS dates derived from the four individual curves are displayed by latitude in Figure 2. According to Figure 2, good consistency can be found among the EOS dates derived from the three indicators except for NDVI. EOS dates derived from NDVI are generally 60 days larger and 40 days larger than those derived from the other three indicators for North America and Eurasia respectively. For SOS, dates derived from the four indicators are much closer. NDVI and APAR give 10-20 days' earlier SOS dates than GPP in both North America and Eurasia. SOS dates based on the NDVI 2 × PAR curves keep consistent with those based on GPP curves across latitudes over years in America. In Eurasia, SOS dates based on the former are generally smaller than those based on the latter, especially at higher latitudes. This is because North America consists mainly of ENF while in Eurasia forests are ENF-DBF mixed. In addition, EOS dates were more strongly dependent on latitude than SOS dates, especially over Eurasian forests. results are expected when using APAR or VI-based LUE model to monitor forests' seasonal dynamics rather than using NDVI itself.
Compared to ENF, DBF shows narrower seasonal cycles (Figure 1g-i). This is due to the fact that in the spring DBF rely on newly grown leaves for photosynthesis while new leaves only give limited contribution to ENF's total photosynthesis, while in the fall DBF experiences more sharp decline of canopy greenness than ENF. Both NDVI and APAR show broader seasonal cycles than GPP. It is shown that for DBF, the NDVI 2 × PAR curves still slightly overestimate the duration of growing season, especially in the spring. This is possibly because photosynthesis for DBF starts only when greenness and temperature increase to a higher level than ENF, leading more lag between spring VIs and GPP. It might be interpreted that it is more difficult to use VIs to track seasonal LUE changes in DBF. One possible explanation lies in that for DBF, seasonal variations of NDVI are jointly affected by changes of canopy structure and leaf pigment while for ENF whose canopy structure shows less seasonal dependence, NDVI variations are mainly explained by changes of canopy pigment concentration, which is directly related to LUE. In this sense, NDVI is more reasonable a proxy for LUE in ENF, where the VI-based LUE model performs better in the present study, than in DBF. Although monitoring phenology over ENF is generally more difficult than over DBF using remote sensing (e.g., [25]), the results in the present study suggest that the VI-based LUE model can track GPP dynamics more precisely in ENF.
SOS and EOS dates derived from the four individual curves are displayed by latitude in Figure  2. According to Figure 2, good consistency can be found among the EOS dates derived from the three indicators except for NDVI. EOS dates derived from NDVI are generally 60 days larger and 40 days larger than those derived from the other three indicators for North America and Eurasia respectively. For SOS, dates derived from the four indicators are much closer. NDVI and APAR give 10-20 days' earlier SOS dates than GPP in both North America and Eurasia. SOS dates based on the NDVI 2 × PAR curves keep consistent with those based on GPP curves across latitudes over years in America. In Eurasia, SOS dates based on the former are generally smaller than those based on the latter, especially at higher latitudes. This is because North America consists mainly of ENF while in Eurasia forests are ENF-DBF mixed. In addition, EOS dates were more strongly dependent on latitude than SOS dates, especially over Eurasian forests.

Discussion
NDVI has been widely used to track seasonal dynamics of photosynthesis and phenology based on NDVI are often compared to both flux tower measurements and remote sensing measurements. Not surprisingly, NDVI-based growing season has been found to be longer than GPP-based growing season. This is because, in the spring, photosynthesis lags behind leaf growth [26] while in the fall, leaf decay lags behind photosynthesis decline [25]. In order to improve VI-based phenology, temperature, which high-latitude forests are sensitive to and therefore dominating their phenology, is often taken into consideration (e.g., [7]). However, for boreal forests whose phenology is mainly determined by environmental factors, seasonal VI cycles also carry information of temperature and precipitation. Instead of introducing environmental stress factors, the present study aims to improve NDVI-based phenology at the perspective of using VI-based GPP estimation model.
Because NDVI has been considered as a proxy for fPAR, the key of using the VI-based LUE model to track seasonal cycles of GPP is whether NDVI can be a reasonable proxy for LUE at seasonal scale. Remote sensing of LUE from greenness-related VIs is challenging due to the fact that LUE is jointly controlled by multiple factors and is sensitive to changes of environmental stress. As a matter of fact, VIs often show obvious lag when environmental stress occurs. At short-term scale (e.g., diurnal), where physiological factors such as non-photochemical quenching (NPQ) dominate LUE changes [15], greenness-related VIs generally fail to track LUE changes [14]. At seasonal scale, however, both LUE and NDVI are sensitive to canopy chlorophyll concentration [27,28]. Thus NDVI may track seasonal changes of LUE. These facts give opportunities for using VI as a proxy for LUE at seasonal scale, especially for phenology study. As a matter of fact, linear relationships between seasonal NDVI and LUE has been reported by empirical analyses [18,19]. Another possible reason for why NDVI-based phenology is improved by applying VI-based LUE model is that NDVI 2 might be a good proxy for ecosystem/apparent LUE (i.e., the ratio of GPP and PAR). Numerically, squaring NDVI gives smaller value and the squared NDVI can be closer to actual fPAR related to GPP due to the fact that only a part of the canopy APAR (PAR absorbed by photosynthetically related pigments) is used for photosynthesis. Since LUE model is conceptual [11], NDVI 2 × PAR might be a direct estimation of GPP regardless of using one NDVI as proxy for fPAR while using another as proxy for LUE.
Results of the present study show that VI-based phenology is improved towards better tracking seasonal photosynthesis (GPP) dynamics through applying a VI-based LUE model. Better performance of the proposed method is found over ENF than DBF. This is because VI and photosynthesis have different mutual influencing mechanisms over these two plant functional types. NDVI variations of ENF are more controlled than those of DBF by changes of canopy pigment concentration, which is the main factor determining canopy LUE. Using remotely sensed VIs to monitor phenology over ENF is generally more difficult because ENF has smaller annual VI fluctuations.
A limitation of this study is that errors of MPI-GPP were not considered although MPI-GPP is generally considered to have overall good quality.

Discussion
NDVI has been widely used to track seasonal dynamics of photosynthesis and phenology based on NDVI are often compared to both flux tower measurements and remote sensing measurements. Not surprisingly, NDVI-based growing season has been found to be longer than GPP-based growing season. This is because, in the spring, photosynthesis lags behind leaf growth [26] while in the fall, leaf decay lags behind photosynthesis decline [25]. In order to improve VI-based phenology, temperature, which high-latitude forests are sensitive to and therefore dominating their phenology, is often taken into consideration (e.g., [7]). However, for boreal forests whose phenology is mainly determined by environmental factors, seasonal VI cycles also carry information of temperature and precipitation. Instead of introducing environmental stress factors, the present study aims to improve NDVI-based phenology at the perspective of using VI-based GPP estimation model.
Because NDVI has been considered as a proxy for fPAR, the key of using the VI-based LUE model to track seasonal cycles of GPP is whether NDVI can be a reasonable proxy for LUE at seasonal scale. Remote sensing of LUE from greenness-related VIs is challenging due to the fact that LUE is jointly controlled by multiple factors and is sensitive to changes of environmental stress. As a matter of fact, VIs often show obvious lag when environmental stress occurs. At short-term scale (e.g., diurnal), where physiological factors such as non-photochemical quenching (NPQ) dominate LUE changes [15], greenness-related VIs generally fail to track LUE changes [14]. At seasonal scale, however, both LUE and NDVI are sensitive to canopy chlorophyll concentration [27,28]. Thus NDVI may track seasonal changes of LUE. These facts give opportunities for using VI as a proxy for LUE at seasonal scale, especially for phenology study. As a matter of fact, linear relationships between seasonal NDVI and LUE has been reported by empirical analyses [18,19]. Another possible reason for why NDVI-based phenology is improved by applying VI-based LUE model is that NDVI 2 might be a good proxy for ecosystem/apparent LUE (i.e., the ratio of GPP and PAR). Numerically, squaring NDVI gives smaller value and the squared NDVI can be closer to actual fPAR related to GPP due to the fact that only a part of the canopy APAR (PAR absorbed by photosynthetically related pigments) is used for photosynthesis. Since LUE model is conceptual [11], NDVI 2 × PAR might be a direct estimation of GPP regardless of using one NDVI as proxy for fPAR while using another as proxy for LUE.
Results of the present study show that VI-based phenology is improved towards better tracking seasonal photosynthesis (GPP) dynamics through applying a VI-based LUE model. Better performance of the proposed method is found over ENF than DBF. This is because VI and photosynthesis have different mutual influencing mechanisms over these two plant functional types. NDVI variations of ENF are more controlled than those of DBF by changes of canopy pigment concentration, which is the main factor determining canopy LUE. Using remotely sensed VIs to monitor phenology over ENF is generally more difficult because ENF has smaller annual VI fluctuations.
A limitation of this study is that errors of MPI-GPP were not considered although MPI-GPP is generally considered to have overall good quality.

Conclusions
In this study, we investigate the feasibility of using the VI-based LUE model to track seasonal dynamics of GPP over high-latitude forests. Our results showed that VI-based phenology of ENF can be significantly improved through applying the VI-based LUE model. One of the potential advantages of using VI-based LUE model to monitor photosynthesis dynamics is that apart from VI, no external data source is required except for PAR, which can be easily obtained from satellite or in-situ observations. Another advantage is that there is no need to estimate absolute LUE or GPP values, neither of which can be directly obtained from remote sensing data. Several recent studies have employed satellite-based solar-induced chlorophyll fluorescence (SIF) to track seasonal photosynthesis cycles, as SIF has been considered as a proxy for GPP [29]. The advantage of using SIF to monitor phenology is that it carries information of PAR, fPAR (greenness) and LUE [30]. However, as satellite retrieval of SIF could be noisy and global long-term SIF products are lacking, VI is still the most widely used indicator for the phenology study. Our results show the potential of using VI-based LUE model to monitor photosynthesis dynamics more precisely.