Phenology-Based Maximum Light Use Efﬁciency for Modeling Gross Primary Production across Typical Terrestrial Ecosystems

: The maximum light use efﬁciency (LUE) ( ε 0 ) is a key essential parameter of the LUE model, and its accurate estimation is crucial for quantifying gross primary production (GPP) and better understanding the global carbon budget. Currently, a comprehensive understanding of the potential of seasonal variations of ε 0 in GPP estimation across different plant functional types (PFTs) is still lacking. In this study, we used a phenology-based strategy for the estimation of ε 0 to ﬁnd the optimal photosynthetic responses of the parameter in different phenological stages. The start and end of growing season (SOS and EOS) from time series vegetation indices and the camera-derived greenness index were extracted across seven PFT ﬂux sites using the methods of the hybrid generalized additive model (HGAM) and double logistic function (DLF). Optimal extractions of SOS and EOS were evaluated, and the ε 0 was estimated from ﬂux site observations during the optimal phenological stages with the light response equation. Coupled with other obligatory parameters of the LUE model, phenology-based GPP (GPP phe-based ) was estimated over 21 site-years and compared with vegetation photosynthesis model (VPM)-based GPP (GPP VPM ) and eddy covariance-measured GPP (GPP EC ). Generally, GPP phe-based basically tracked both the seasonal dynamics and inter-annual variation of GPP EC well, especially at forest, cropland, and wetland ﬂux sites. The R 2 between GPP phe-based and GPP EC was stable between 0.85 and 0.95 in forest ecosystems, between 0.75 and 0.85 in cropland ecosystems, and around 0.9 in wetland ecosystems. Furthermore, we found that GPP phe-based was signiﬁcantly improved compared to GPP VPM in cropland, grassland, and wetland ecosystems, implying that phenology-based ε 0 is more appropriate in the GPP estimation of herbaceous plants. In addition, we found that GPP phe-based was signiﬁcantly improved over GPP VPM in cropland, grassland, and wetland ecosystems, and the R 2 between GPP phe-based and GPP EC was improved by up to 0.11 in cropland ecosystems and 0.05 in wetland ecosystems compared to GPP VPM , and RMSE was reduced by up to 5.90 and 2.11 g C m − 2 8 day − 1 , respectively, implying that phenology-based ε 0 in herbaceous plants is more appropriate for GPP estimation. This work highlights the potential of phenology-based ε 0 in understanding the seasonal variation of vegetation photosynthesis and production.


Introduction
Gross primary productivity (GPP) of the terrestrial ecosystem, which reflects the carbon dioxide (CO 2 ) flux fixed by plants through the process of photosynthesis, is a crucial factor in the global carbon cycle. Conventionally, the eddy covariance (EC) technique is recognized as one of the best micrometeorological methods for measuring the net exchange (NEE) of CO 2 between the atmosphere and the surface of various ecosystems and ecosystem respiration (ER) across diverse terrestrial ecosystems. The observed NEE, partitioned into GPP and ecosystem respiration, can be used for indirect GPP estimation [1]. EC flux measurements from continental or global flux networks provide fundamental data for calibrating and validating the parameters of ecological models to improve GPP estimation capability, facilitating numerous studies of carbon fluxes over various terrestrial ecosystems [2]. These flux sites' data integrated with satellite remote sensing and climate data have been widely used to support the development of GPP estimation models at the site, regional, and global scales [3,4].
Light use efficiency (LUE) models, driven by vegetation indices (e.g., the enhanced vegetation index, EVI) obtained from remotely sensed imagery and climate data, have been extensively used for estimating GPP from site scale to global scale due to their advantages of a well-accepted explanation, few input parameters, and ease of operation. These models estimate GPP as a product of absorbed photosynthetically active radiation (APAR) and LUE (ε g ) (GPP = APAR × ε g ), which is derived from downscaling the maximum LUE (LUE max , denoted as ε 0 ) by the scaling of environmental constraints (e.g., temperature, water, CO 2 , etc.). Generally, the ε 0 is initially regarded as a constant for global vegetation in LUE models, such as C-Fix [5] and EC-LUE [6], and then considered as a biome-specific constant, as in MOD17 [7], or a constant dependent on product types of carbon (namely C3 and C4), as in VPM and TEC models [3,8], as well as a two-leaf-specific constant in the TL-LUE model [9]. However, growing ecological communities think of ε 0 as a dynamic value rather than a constant [10,11].
In reality, ε 0 was significantly influenced by vegetation canopy absorption of solar radiation [12] and the internal physiological characteristics (e.g., chlorophyll content) of plants [13] in earlier studies. Based on the fact that there is an essential increase in light use efficiency under overcast conditions, a cloudiness index was used to dynamically regulate ε 0 in CFLUX and CI-LUE models [12,14]. Considering the nonlinear response of vegetation photosynthesis to solar radiation induced by vegetation photosynthesis varies with the dynamics of solar radiation, Xie et al. (2023) [10] proposed a PAR-regulated dynamic ε 0 in GPP estimation based on the LUE model, while several recent studies have shown that ε 0 greatly varies during the vegetation growing season [15] as a result of the varying correlation between chlorophyll content and photosynthesis efficiency at various growing stages [16]. Dynamic ε 0 was parameterized in the phenological transitions of paddy rice [17] and the leaf area index (LAI)-based green-up stage and senescence stage [11] at flux sites to improve agroecosystem GPP estimation. So far, studies on the comprehensive performances of phenology-based dynamic ε 0 for GPP estimation in broader terrestrial ecosystems (such as forests, grasslands, and wetlands) are few, but performance evaluation can provide a deep understanding of dynamic ε 0 in GPP estimation.
More specifically, Huang et al. (2021) [17] estimated ε 0 during four phenological stages corresponding to the physiological features of paddy rice. The results indicated that phenology-based ε 0 is more appropriate for GPP estimation in the paddy field ecosystem. In another relevant study, Huang et al. (2022) [11] found seasonal variations of ε 0 in agroecosystems and re-parameterized ε 0 in the green-up and senescence stages to provide high-accuracy GPP estimates. These studies have demonstrated that phenology-based ε 0 performs better than static or constant ε 0 in GPP estimation because it is more consistent with phenology variations in the vegetation. In terms of the global terrestrial ecosystem, phenological stages diverge as a result of the changes in photosynthesis efficiency and rate caused by plant physiological characteristics [10,11], and they are appropriated for estimating ε 0 in vegetation phenological stages.
Remote Sens. 2023, 15, 4002 3 of 23 In order to evaluate the improvement of phenology-based ε 0 in GPP estimates across different terrestrial ecosystems, the extraction of identical phenology indicators is a key step. The phenological stage indices, including the start (SOS), end (EOS), and length (LOS) of the growing season, have been successfully derived from remotely sensed imagery or ground observations with various algorithms [18]. Recent studies have showed that extraction algorithms such as the hybrid generalized additive model (HGAM) and double logistic function (DLF) are robust and appropriate for estimating SOS/EOS/LOS from satellite normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and in situ GPP observations [19][20][21]. Moreover, automated digital camera imagery collected through the Phenology Camera (PhenoCam) observation network can provide 'near-surface' observations of vegetation phenology with high temporal resolution [22,23]. These PhenoCam images have been successfully used to evaluate satellite-derived phenology [24][25][26] and test phenology algorithms in diverse ecosystems [27][28][29]. Therefore, a thorough evaluation of the implications of phenology extraction algorithms on phenology-based ε 0 and a comprehensive comparison of phenology-based ε 0 across different ecosystems are required.
In an effort to address the issue of the current ε 0 estimation approach not fully taking into account the identical phenological stages of various PFTs, we proposed a framework for phenology-based regulation of ε 0 and made a comparison of GPP estimates from phenologybased ε 0 . Specifically, the objectives of the study were to: (1) compare phenology-based ε 0 estimated from NDVI, EVI, and camera-derived greenness index based on two well-known algorithms (HGAM and DLF) among forest, grassland, cropland, and wetland ecosystems and (2) assess the degree of improvement in GPP estimation by use of phenology-based ε 0 in the LUE model.

FLUXNET Data
In this study, data from seven flux sites from the FLUXNET dataset were collected, covering evergreen needleleaf forest, deciduous broadleaved forest, paddy field, dryland, grassland, and wetland ecosystems. The criterion used to select flux sites was that they had to contain observations that include both flux data and phenology camera records from the same time period ('Data Availability' in Table 1). Therefore, we selected the time scale of these flux tower sites as 2010 to 2015. The FLUXNET provides two public datasets: the FLUXNET2015 Dataset and the FLUXNET-CH4 Community Product. To improve consistency and comparability, both datasets offer preprocessing and qualitycontrolled flux data using consistent methodologies [30][31][32]. The variables used in this study included gross primary productivity ('GPP_NT_VUT_REF'), net ecosystem exchange ('NEE_VUT_REF'), ecosystem respiration ('RECO_VUT_REF'), air temperature ('TA_F'), and photosynthetic photon flux density incident ('PPFD_IN'). The variables of the USTAR ('VUT') threshold for selected years were used in this study [33][34][35], and detailed data processing method can be found on the FLUXNET website (https://fluxnet.org/data/ fluxnet2015-dataset/data-processing (accessed on 22 July 2023)). High-quality data were refined to remove abnormal data with PPFD_IN_QC values of −9999 or 0. With regard to GPP, the variables of 'GPP_NT_VUT_REF' were recognized as reference daily GPP (g C·m −2 ·d −1 ). The variables of half-hour and hourly input were scaled to a daily scale for calibration and validation of the model estimates.
Howland Research Forest is located in Piscataquis County, central Maine. It is characterized by evergreen coniferous forest, with soils formed on coarse-grained granite. The vegetation consists of approximately 66% red spruce, 25% eastern hemlock, and 34% other coniferous trees and hardwoods [36]. Harvard Forest is situated in Petersham County, Massachusetts. It is a deciduous broadleaf forest, with soils consisting of rocky sandy loam and glacial till. The dominant vegetation types include red oak, sugar maple, eastern white pine, and yellow birch [21]. Morgan Monroe State Forest is located in Morgan and Monroe Counties, Indiana. It is a deciduous broadleaf forest characterized by a low-relief landscape with ridges and ravines, where the maximum relief is less than 60 m. The main vegetation types are sugar maple, tulip poplar, yellow poplar, and various oak species [37]. US-Ne2 is situated near Mead, Nebraska, at the University of Nebraska Agricultural Research and Development Center. It represents a typical corn-soybean rotation agricultural ecosystem, with corn planted in odd-numbered years. Unlike US-Ne3, which relies solely on rainfall, US-Ne2 is equipped with a central-pivot irrigation system [38]. Twitchell Island rice field is located within the Sacramento-San Joaquin Delta, near the Pacific coast. This rice field ecosystem experiences a Mediterranean climate with distinct wet and hot seasons [39]. Mayberry Wetland is situated on Sherman Island and is a constructed restored wetland built using the heterogeneity trenching method. The wetland vegetation is predominantly cattail, and it exhibits the highest heterogeneity in water depth among the restored wetlands from the same period [40,41]. Vaira Ranch is located near Ione, California. The primary vegetation consists of annual C3 grassland, influenced by a Mediterranean climate. These annual grasslands exhibit higher activity during the warm and moist winter season compared to during the summer season [42]. The distribution of the above sites is shown in Figure 1. Site US-Ne2 is covered by dryland, and site US-Twt is covered by paddy field. T min , T max and T opt are the minimum temperature, the maximum temperature and the optimal temperature, which were used in previous studies based on site observations.

Phenology Camera Data
The PhenoCam tracks vegetation phenology in various ecosystems, initially in North America and then globally [43,44]. Cameras are pointed in the dominant wind direction in order to match the eddy covariance tower footprint [45]. Digital images are stored as JPG files (image resolution of about 1900 × 960 pixels, three color channels with 8-bit RGB color information) at half-hourly intervals for several hours per day with exposure mode. If digital images contain devices or the sky, several regions of interest (ROIs) are used to extract vegetation enclosing all areas in the foreground. As GCC (green chromatic coordinate) is widely used for monitoring vegetation canopy development and quantifying canopy greenness [45,46], it is estimated from the following equation: where the subscripts DN G , DN R , and DN B represent the red, green, and blue channels, respectively. The daily GCC is smoothed by a filter based on the official algorithm to reduce some uncertainty [47] and aggregated at 8-day intervals to match EVI temporal resolution for extracting phenology. Because there were no time periods overlapping the FLUXNET and PhenoCam databases, the GCC from the US-Ne2 site (dryland environment) could not be used in this analysis.
stored wetland built using the heterogeneity trenching method. The wetland vegetation is predominantly cattail, and it exhibits the highest heterogeneity in water depth among the restored wetlands from the same period [40,41]. Vaira Ranch is located near Ione, California. The primary vegetation consists of annual C3 grassland, influenced by a Mediterranean climate. These annual grasslands exhibit higher activity during the warm and moist winter season compared to during the summer season [42]. The distribution of the above sites is shown in Figure 1.

MODIS Data
The MODIS surface spectral reflectance (MOD09A1 V6), FPAR, and leaf area index (LAI) products (MCD15A3H V6) are accessible at a spatial resolution of 500 m on the Google Earth Engine (GEE) platform and were used in the study. The MOD09A1 and MCD15A3H products are 8-day and 4-day composite data which choose the best pixels filtered by their quality labels. We used red (620-670 nm), blue (459-479 nm), near infrared (NIR, 841-876 nm), and shortwave infrared (SWIR, 2105-2155 nm) bands to calculate NDVI, EVI, and LSWI (land surface water index) [48][49][50] in the following equations: where ρ Blue , ρ Red , ρ NIR , and ρ SWIR are the spectral reflectance of the blue, red, NIR, and SWIR bands, respectively. G is set to 2.5, C 1 is set to 6, C 2 is set to 7.5, and L is set to 1 in Equation (3). The MODIS sensor may have a small amount of abnormal data or missing data due to its internal components or external solar radiation, whose time series data need to be gap-filled and smoothed through interpolation and filter. The Savitzky-Golay (SG) filter employs a locally weighted regression to fit the data within the window as a polynomial. It replaces the original value at the center point of the window with the polynomial value. The degree of fitting of the SG filter depends on the window size and the order of the polynomial. The sliding window sizes for SG filtering of NDVI, EVI, LSWI, and LAI time series data were set to 11 (for a few sites where data anomalies still existed after processing, the window sizes were set to 15), and the order of the polynomial was set to 3. The 4-day LAI products were averaged every two observations to match the time resolution of MOD09A1. Likewise, we filtered daily GCC data and took the maximum value of GCC within every 8-day period as the value for that period.

Methodology
The overall framework of this study is shown in Figure 2. Some specific data processing methods and GPP estimation models are detailed and described in this section.
window with the polynomial value. The degree of fitting of the SG filter depends on the window size and the order of the polynomial. The sliding window sizes for SG filtering of NDVI, EVI, LSWI, and LAI time series data were set to 11 (for a few sites where data anomalies still existed after processing, the window sizes were set to 15), and the order of the polynomial was set to 3. The 4-day LAI products were averaged every two observations to match the time resolution of MOD09A1. Likewise, we filtered daily GCC data and took the maximum value of GCC within every 8-day period as the value for that period.

Methodology
The overall framework of this study is shown in Figure 2. Some specific data processing methods and GPP estimation models are detailed and described in this section.

Figure 2.
The overall framework of the study.

Double Logistic Function
The DLF method has been fast developed and widely used to extract phenological indices at site or regional scale. It is a common nonlinear function that is made up of two logistic functions, each with an inflection point and a curvature parameter. We used the following seven-parameter logistic function to extract SOS and EOS at the seven flux sites. The parameters of the DLF can be adjusted to fit various ecosystem responses, making it highly suitable for different environmental conditions and ecosystem types [19,51].

Double Logistic Function
The DLF method has been fast developed and widely used to extract phenological indices at site or regional scale. It is a common nonlinear function that is made up of two logistic functions, each with an inflection point and a curvature parameter. We used the following seven-parameter logistic function to extract SOS and EOS at the seven flux sites. The parameters of the DLF can be adjusted to fit various ecosystem responses, making it highly suitable for different environmental conditions and ecosystem types [19,51].
where Y(t) represents the observed data value at time t on DOY (day of year). α 1 is the value during the winter dormant period, α 2 -α 1 is the amplitude between the winter dormant period and the spring and early-summer peak, and α 3 -α 1 is the amplitude between the winter dormant period and the late summer and autumn peak. ∂ 1 and ∂ 2 are curvature metrics to adjust the slope of the logistic function. β 1 and β 2 are the inflection points in DOY of the transitions for vegetation green-up and senescence stages. In this study, the fitted parameters of β 1 and β 2 represent the SOS and EOS, and the length of the growing season is the difference between β 1 and β 2 .

Hybrid Generalized Additive Model
The HGAM is a statistical model that combines the advantages of the generalized additive model (GAM) and hybrid model which was proposed for estimating SOS and PPY (peak photosynthesis timing) from GPP data by Yang et al. [20]. The GAM is a nonlinear modeling approach that describes the nonlinear effects of independent variables as smooth functions. Taking into account the ambiguous representation of long time series feature points from raw discrete data, it was developed into a hybrid model. A detailed description of the HGAM can be found in Yang et al. [20]. When modeling NDVI, EVI, and GCC time series data with the GammaGAM class in the PyGAM library, the data need to be preprocessed with gap filling and smoothing to make the input data more stable and acceptable after entering the model. The cubic spline interpolation was used for gap filling, and the smoothing operation is explained in Section 2.3.
This study employed ten splines to model the relationship between independent and response variables. Specifically, each spline function defined a spline node, which divided the range of independent variable into ten intervals. The interval ranges were divergent, making the modeling procedure more flexible. Compared with the DLF method, the HGAM does not have a concept similar to the inflection point and uses local tuning thresholds to determine the temporal indices of the key phenological stages. Since the range of GCC derived from the camera images of the flux sites was usually less than 0.2 (mostly between 0.3 and 0.5), the mean of the maximum and minimum values (the difference between the maximum and minimum values of GCC data for each group of experiments in the study was less than 0.2) was used as the tuning threshold rather than the averaged amplitude. Then, we set the threshold in the smoothed NDVI, EVI, and GCC curve to identify the SOS and EOS. Based on the SOS and EOS estimated from the vegetation indices, three phenological stages (i.e., before SOS, SOS-EOS, and after EOS) were derived for GPP estimation in the next step.

Structure of the LUE Model
The general structure of a light use efficiency model to estimate GPP can be formulated as follows: where ε g is the actual maximum light use efficiency (µmol CO 2 /µmol photosynthetic photon flux density, PPFD), PAR is the incident photosynthetically active radiation (µmol PPFD), and FAPAR is the fraction of PAR absorbed by vegetation. The ε g is affected by temperature and water stress and can be expressed as: where ε 0 is the apparent quantum yield or maximum LUE (µmol CO 2 /µmol PPFD or g C/mol PPFD), and T scalar and W scalar are the scalars for the effects of temperature and water on the light use efficiency of vegetation, respectively. Using the equation developed for the terrestrial ecosystem model [52], T scalar from the vegetation photosynthesis model (VPM) was used in this study [53]: where T min , T max , and T opt denote the minimum, maximum, and optimal temperature for photosynthetic activity, respectively. Their values were determined based on the previous research on GPP estimation in specific ecosystems, along with the sites' real circumstances [3]. W scalar was estimated from a satellite-derived, water-sensitive vegetation index: (12) where LSWI max is the maximum land surface water index (LSWI) during the vegetation growing seasons and was calculated separately for each year. Generally, ε 0 is dependent on the PFTs and can be estimated from light-response curves to fit time series data of net ecosystem exchange (NEE), photosynthetically active radiation (PAR), and ecosystem respiration (Re). Specifically, this study employed a rectangular hyperbolic function (the Michaelis-Menten light response equation) to fit the three parameters mentioned above from eddy covariance measurements [54]: where GEE max is the maximum rate of ecosystem gross photosynthesis (µmol CO 2 ·m −2 ·s −1 ), and Re is the daytime ecosystem respiration. As ε 0 values vary with phenological stages of different PFTs [17,55], we fitted ε 0 during three phenological stages (i.e., before SOS, SOS-EOS, and after EOS), which were extracted from the time series vegetation indices based on the algorithms of DLF and HGAM (see detailed description in Section 3.1).

The Calculation of FAPAR
Some earlier research on LUE models employed the fraction of PAR absorbed by the vegetation canopy (FAPAR canopy ) to estimate APAR, using a linear relationship between FAPAR canopy and VIs as well as LAI as FAPAR canopy proxies [5,8,56]. Subsequently, as communities discovered that the fraction of PAR absorbed by chlorophyll (FAPAR chl ) or photosynthetically active vegetation (PAV, FAPAR pav ) contributes more to vegetation photosynthesis, FPAR chl or FAPAR pav were found to be appropriate for estimating GPP. In this study, we selected four common relationships between FAPAR chl /FAPAR canopy and satellite-based vegetation indicators as follows: FAPAR chl3 = min(max(d × NDVI + e, 0), 1) FAPAR canopy = 1 − e k×LAI (17) where Equations (14) and (15) are linear functions between EVI and FAPAR chl where a is set to 1 [53], and b and c are, respectively, 0.1 and 1.25 [10]. Equation (16) removes the influence of bare soil on NDVI, where d and e are set to 1.24 and −0.168 [57]. k is the extinction coefficient set to 0.5 in Equation (17) [58].

Evaluation of Phenological Stages and GPP Estimation
In order to evaluate the performance of the DLF and HGAM in phenology stages extraction, we compared modeled SOS and EOS with the digital images acquired on the same day, as well as 8 days ahead of and lagging behind the date, by visual inspection. The phenology cameras were able to capture the seasonal variations in foliage pigments of most PFTs, which result in an increase in canopy greenness in the spring and a decrease in canopy greenness in the fall. The agreement between the 8-day estimated GPP and observed GPP was evaluated by the coefficient of determination (R 2 ) and the root-mean-square error (RMSE). Furthermore, GPP estimates from the VPM model were used as a reference to evaluate phenology-based modeled GPP. ε 0 and T scalar were set according to site-specific references [3,6,39,53], and others were default parameters.
In the equation, n represents the total number of samples, Y pre_t is the predicted value for the current period, Y obs_t is the observed value for the current period, Y pre_mean is the mean value of the predicted values, and Y obs_mean is the mean value of the observed values.

SOS and EOS Estimated from Vegetation Indices Based on DLF and HGAM at Flux Sites and Their Validation
The SOS and EOS derived from the DLF and HGAM models varied with the time series observations. Among the 60 experiments, the HGAM method successfully extracted the SOS and EOS within the effective time range (referring to 46 eight-day periods in a year) in all cases. However, the DLF method only extracted the SOS and EOS within the effective time range in 45 out of the 60 experiments. The calculated effective extraction rate for the HGAM was 100%, while, for DLF, it was 75%. Notably, among the 18 experiments using GCC, DLF failed to provide valid results in 12 cases. Therefore, despite the different vegetation indices used in phenology extraction of different ecosystems, the HGAM consistently produced valid results and exhibited better robustness compared to DLF.
Specific dates of SOS and EOS in the DOY with an 8-day interval (the same as used below) are summarized in Table 2 We used phenology camera photos to evaluate the accuracy of the key phenological events. The photos acquired at the SOS and EOS extracted from NDVI, EVI, and GCC using the DLF and HGAM models were downloaded from PhenoCam. A red rectangular box enclosing the vegetation area in the foreground was set in each photo for visual inspection. We took the visual interpretation of US-MMS (in 2014) and US-Myb (in 2013) as examples to represent the woody plant and herbaceous plant. At US-MMS, green leaves were less visible in the photos where the acquired date of SOS was estimated from NDVI_DLF, NDVI_HGAM, EVI_HGAM, and GCC_HGAM, while many new green leaves had emerged in the photo whose acquired date of SOS was estimated from EVI_HGAM. The EOS date derived from the NDVI_HGAM, GCC_HGAM, and NDVI_DLF was obviously lagging behind, which could be seen by the large number of yellow leaves or chlorophyll-degraded leaves in these photos (Figure 3a). The photo witnessed the forest growing stage when the camera date was derived from the SOS of the EVI_DLF, and the EOS estimated from the EVI_DLF and EVI_HGAM matched the time that the vegetation started to turn yellow in the camera images well. Overall, the SOS and EOS obtained by EVI_HGAM were more compatible with the vegetation growth stages, and their dates were very close to those in the research conclusions of Cheng et al. [59]. With regards to the restored wetland of US-Myb in 2013, the SOS was estimated to be in the 12th 8-day interval from EVI_DLF, NDVI_HGAM, and EVI_HGAM, while it was estimated from GCC_HGAM to be in the 15th 8-day interval. As shown in Figure 3b, the vegetation was still yellow in the 12th 8-day interval and started to green up in the 15th 8-day interval, which basically coincided with the SOS in DOY 105-140 in Fang et al. [60]. Similar evaluations of estimates of SOS and EOS were implemented at other sites. Specifically, the EVI_HGAM was the appropriate model to extract SOS and EOS at US-Ho1, US-Ha1, US-MMS, US-Ne2, and US-Twt, GCC_HGAM was suitable for US-Myb, and EVI_DLF performed better than other models in US-Var (Supplementary Materials, Figure S2.1-S2.21).  (16.25) vegetation was still yellow in the 12th 8-day interval and started to green up in the 15th 8day interval, which basically coincided with the SOS in DOY 105-140 in Fang et al. [60]. Similar evaluations of estimates of SOS and EOS were implemented at other sites. Specifically, the EVI_HGAM was the appropriate model to extract SOS and EOS at US-Ho1, US-Ha1, US-MMS, US-Ne2, and US-Twt, GCC_HGAM was suitable for US-Myb, and EVI_DLF performed better than other models in US-Var (Supplementary Materials, Figure  S2

Estimates of the Maximum LUE (ε 0 ) during Phenological Stages
The phenology-based maximum light use efficiency (ε 0 ) at the seven flux sites is shown in Table 3. Phenological stage 1 (PS1) refers to the period from the first day in the DOY to the SOS, phenological stage 2 (PS2) refers to the second period from the SOS to the EOS, and phenological stage 3 (PS3) refers to the third period from the EOS to the last day in the DOY. At the US-Ho1 site from 2013 to 2015, the maximum ε 0 appeared in the PS3 and varied from 0.042 to 0.049 µmol CO 2 /µmol PPFD, while the minimum ε 0 appeared in the PS1 and differed from 0.030 to 0.036 µmol CO 2 /µmol PPFD. For the DBF site, ε 0 showed considerable ranges from 0.037 to 0.062 µmol CO 2 /µmol PPFD at US-Ha1 and 0.013-0.040 µmol CO 2 /µmol PPFD at US-MMS. In regard to the cropland site of US-Ne2 and US-Twt, the maximum ε 0 was synchronous in the PS2 with the value of 0.065 in the dryland and 0.045 in the paddy field. The ε 0 of the wetland (US-Myb site) in PS2 was relatively higher than that in PS1 and PS3, with the maximum value being approximate to 0.38 from 2012 to 2014. For the grassland site of US-Var during 2012-2014, the ε 0 showed substantial variations among the PS1, PS2, and PS3, especially the minimum value of 0.01 in PS1 and the maximum value of 0.084 in PS2 in the year of 2014. The maximum ε 0 ranged from 0.042 µmol CO 2 /µmol PPFD to 0.084 µmol CO 2 /µmol PPFD in the PS2 of the three years.

Comparison of GPP Estimates with Phenology-Regulated ε 0 and Different FAPAR Proxies
A total of 21 site-years GPP values were simulated with phenology-based ε 0 and four FAPAR proxies. For one site-year GPP simulation, the specific parameters used to generate five GPP estimates are shown in Table 4. GPP5 and GPP EC represent VPM GPP and observed GPP estimated from each flux site. Table 3. Specific values of ε 0 (µmol CO 2 /µmol PPFD) estimated from phenological stages at the seven sites.

Site_ID
Year The temporal dynamics of five GPP products are plotted against the temporal dynamics of GPP EC in Figure 4. At US-Ho1 and US-Ha1, GPP2, GPP3, and GPP5 accurately tracked both the seasonal dynamics and inter-annual variation of GPP EC , while GPP1 and GPP4 somewhat to moderately overestimated GPP EC during these years. At US-MMS, five GPP products overestimated the observed GPP from 2012; however, in 2013 and 2014, GPP2 and GPP3 traced the variability of GPP EC on seasonal and inter-annual scales well. At US-Ne2, five GPP products followed GPP EC in 2010 and 2011 well, whereas these modeled GPP values slightly lagged behind the GPP EC in 2012. The same situation occurred at US-Twt in 2013. In 2012 and 2014, GPP4 slightly overestimated GPP EC . GPP1, GPP2, GPP3, and GPP5 overall fitted the observed GPP in 2012, and GPP2, GPP3, and GPP5 somewhat underestimated GPP EC in 2014. From 2012 to 2014, GPP1 and GPP4 were greatly larger than GPPEC at US-Myb. GPP2, GPP3, and GPP5 generally properly matched with GPP EC in 2013 and 2014, while the three GPP products were moderately smaller than the observed GPP in 2012. At US-Var, GPP1 and GPP4 significantly overestimated GPP EC from 2012 to 2014. In particular, in 2013, the five GPP products moderately lagged behind the observed GPP. and GPP5 overall fitted the observed GPP in 2012, and GPP2, GPP3, and GPP5 somewhat underestimated GPPEC in 2014. From 2012 to 2014, GPP1 and GPP4 were greatly larger than GPPEC at US-Myb. GPP2, GPP3, and GPP5 generally properly matched with GPPEC in 2013 and 2014, while the three GPP products were moderately smaller than the observed GPP in 2012. At US-Var, GPP1 and GPP4 significantly overestimated GPPEC from 2012 to 2014. In particular, in 2013, the five GPP products moderately lagged behind the observed GPP.

SOS/EOS Estimating Methods of Different PFTs
Extracting the phenology of different PFTs from time series satellite data is a preliminary step in ecological remote sensing research. In this study, we used two methods, the HGAM and DLF, on three vegetation indices (EVI/NDVI/GCC) to extract the crucial phenological stages (SOS/EOS) of different vegetation types [61,62].
On the one hand, the three vegetation indices have inherent differences in their characteristics. NDVI serves as the foundation of remote sensing phenology and is sensitive to the biochemical and physiological characteristics of vegetation, commonly used for distinguishing between vegetation and bare soil [48]. Using the dynamic threshold method to extract phenological parameters from NDVI time series can improve the prediction of soil organic carbon content during crop rotation [63]. EVI, as an enhanced vegetation index, exhibits higher sensitivity in areas with high biomass. Research has shown that EVI outperforms NDVI in vegetation monitoring, particularly in regions dominated by evergreen broadleaf forests [49]. Using the sigmoidal logistic function, the mutual comparison and interpretation of phenological indicators derived from MODIS EVI and flux tower GPP were validated in southern Taihang Mountains, China [64]. GCC, a derivative product that was developed after the emergence of phenocamera data, holds great potential

SOS/EOS Estimating Methods of Different PFTs
Extracting the phenology of different PFTs from time series satellite data is a preliminary step in ecological remote sensing research. In this study, we used two methods, the HGAM and DLF, on three vegetation indices (EVI/NDVI/GCC) to extract the crucial phenological stages (SOS/EOS) of different vegetation types [61,62].
On the one hand, the three vegetation indices have inherent differences in their characteristics. NDVI serves as the foundation of remote sensing phenology and is sensitive to the biochemical and physiological characteristics of vegetation, commonly used for distinguishing between vegetation and bare soil [48]. Using the dynamic threshold method to extract phenological parameters from NDVI time series can improve the prediction of soil organic carbon content during crop rotation [63]. EVI, as an enhanced vegetation index, exhibits higher sensitivity in areas with high biomass. Research has shown that EVI outperforms NDVI in vegetation monitoring, particularly in regions dominated by evergreen broadleaf forests [49]. Using the sigmoidal logistic function, the mutual comparison and interpretation of phenological indicators derived from MODIS EVI and flux tower GPP were validated in southern Taihang Mountains, China [64]. GCC, a derivative product that was developed after the emergence of phenocamera data, holds great potential for detecting vegetation canopy development [65]. The study conducted by Fan et al. provided a comprehensive perspective on the phenology of complex wetland landscapes by utilizing the DLF method and GCC data [60]. The amplitude of NDVI variation within a year is frequently greater than that of EVI and GCC due to its sensitive response to the change of vegetation canopy compared to other indices. For the estimation of the SOS and EOS of the forest ecosystem in our study, NDVI-based DLF/HGAM algorithms were inferior to the same EVI-based algorithms.
On the other hand, the principles of the two phenology extraction algorithms are different. The double logistic function algorithm evolved from the simple logistic (SL) method with the inclusion of additional variables that are assumed to better detect the changes in greenness from extreme events [51]. The DLF algorithm has shown its potential in monitoring vegetation seasonal growth dynamics as it showed good performance in fitting growth curves in some PFTs (e.g., it fitted well in ENF, paddy fields, and wetlands in our study). However, the maximum rate of change (both increasing and decreasing) of the fitted VI curves supposed to be the SOS and EOS was not accepted. Adversely, over-fitting in green-up or senescence stages, resulting in false inflection points, may have been recognized as the SOS and EOS. Moreover, we found under-fitting of the EVI curve in growth peaks at US-MMS (DBF site), as well as the identical circumstance for the NDVI curves at US-Ne2 (cropland site) and US-Ho1 (ENF site), which implies that there are large uncertainties relating to DLF algorithms applied at regional scale. Generally, the HGAM algorithms showed better performance than DLF in the extraction of the SOS and EOS of different PFTs. It could capture not only a single peak in cropland/grassland/ENF ecosystems, but also asymmetrical peaks in DBF ecosystems during the peak season. Furthermore, VI curves in the periods before green up and after senescence were well tracked by the HGAM algorithms. The HGAM approach allows for the automatic derivation of the predictor functions, which eliminates the requirement to know the fitting parameters and the ultimate kind of predictive functions in advance [20,61]. This can detect intra-annual variations of NDVI/EVI/GCC in the entire growth season across different PFTs using regression splines and smoothing splines [20]. Phenology camera images, which benefit from daily acquisition and in situ near-surface observation, were more welcome in recent studies of GPP estimation. Our findings demonstrated that a GCC time series derived from phenology camera images was capable of extracting phenological metrics with the HGAM, especially for the wetland ecosystem. This conclusion agrees with Knox et al.'s conclusion [45]: that camera-based indices have significant promise in capturing wetland seasonal fluctuations. It should be noted that the estimation of camera-based indices and visual inspection from camera photos can be more ambiguous given the varying growth status of the vegetation in the foreground and background.

Phenology-Based ε 0 of Different PFTs
During various phenological stages, the photosynthetic responses of different PFTs vary quickly or slowly in response to the climate or environment changes. In recent studies, the ε 0 during four phenological stages of paddy rice was investigated and showed 0.020-0.065 mol CO 2 /mol PPFD in Huang et al. [17], 2021, and 0.014-0.0397 mol CO 2 /mol PPFD in Cheng et al. [55] (an approximate conversion of 2.1 between MJ and mol PPFD was used) with different parameter calibration methods. While we used the same calibration method as Huang et al. to estimate ε 0 , and its values ranged from 0.004 to 0.045 mol CO 2 /mol PPFD during the three phenological stages (before SOS, SOS-EOS, after EOS), we concluded that the ε 0 of rice was comparable even if it was derived from different phenology extraction methods (SL in Huang et al. vs. HGAM in this study). Another new piece of research on the ε 0 of dry crops illustrated that ε 0 starts to swiftly increase during the green-up stage and gradually drops during the senescence stage [11], which is consistent with our finding that there is synchrony between ε 0 and crop growing status. Zhu et al. [54] reported that the ε 0 varied with different grassland types, ranging from 0.029 to 0.102 mol CO 2 /mol PPFD across temperate/tropical/alpine grasslands. In that study, the ε 0 of US-Var estimated from M-M equations between 2001 and 2014 was~0.03 mol CO 2 /mol PPFD, which is analogous to our ε 0 estimates at US-Var between 2012 and 2014. Comparing the dynamic range of ε 0 (0.01-0.084 mol CO 2 /mol PPFD) with a fixed constant ε 0 can remind us to pay more attention when using them in grassland ecosystem. Overall, PFTs have a large annual seasonal variation in canopy structure [66], resulting in dynamic ε 0 during the growing season. The estimation of ε 0 with a phenology-based strategy was proposed to find the optimal response of the parameter to photosynthesis in different phenological periods.

Phenology-Based Methods in GPP Estimation with LUE Models-Strengths and Limitations
Comparisons of four GPP products based on phenology-based ε 0 (denoted as GPP phe-based ) and VPM GPP (GPP VPM ) with GPP EC are shown in Figure 5. For the majority of the flux sites, the overall performances revealed good agreements between the four GPP products and GPP EC , and the variation in accuracy was dependent on the use of FAPAR proxies. For forest ecosystems, FAPAR data products based on NDVI and EVI showed weaker improvement than LAI-based FAPAR in most years. In addition, compared with GPP EC , these GPP phe-based values were comparable to GPP VPM at the forest sites, which implies that the use of phenology-based ε 0 in LUE models has limited effect on improving the accuracy of GPP estimation. For herbaceous plants, the correlations between GPP phe-based and GPP EC were subject to different PFTs. Generally, EVI-and NDVI-based FAPAR was found to be more effective than LAI-based FAPAR in GPP phe-based estimation for grassland and wetland ecosystems. The optimal GPP phe-based estimated from phenology-based ε 0 and EVI/NDVI/LAI-based FAPAR performed better than GPP VPM (0.62 vs. 0.14 in R 2 max), indicating significant enhancement by the phenology-based strategy in the estimation of ε 0 and its advantage of reduced uncertainty of the parameter. In terms of the cropland ecosystem, the LAI-based FAPAR used for GPP phe-based estimation worked better than the others in most years (four out of six site-years). The GPP phe-based was more accurate than the GPP VPM in the rice ecosystem (US-Twt site), but there was not much of difference between them in the dryland ecosystem (US-Ne2 site). Considering the positive influences of dynamic ε 0 based on phenological stages in GPP estimation in paddy field environments from previous research and this work [17,55], ε 0 refined at fine spatial scales is a practicable way to extrapolate site GPP to regional or global scale for rice ecosystems.
We compared vegetation photosynthesis model (VPM) GPP (GPP VPM ), phenologybased VPM GPP (GPP phe-based ), and eddy covariance-measured GPP (GPP EC ) across five ecosystems. The substantial improvement in GPP reflects that ε 0 is influenced by the vegetation canopy absorption of solar radiation and internal physiological traits of plants. Under different solar radiation conditions, ε 0 showed different photosynthesis efficiency and rates, for example, increases in light use efficiency under overcast conditions and saturation of vegetation photosynthesis to solar radiation under high solar radiation. Additional factors to be consider are variations in the physiological traits of the vegetation itself, such as changes in chlorophyll content. The distinctive tillage characteristics of drylands and paddy fields, the instability of artificially restored wetlands, and the unique climatic circumstances of grasslands bring about significant changes in these ecosystems compared to in forest ecosystems. Therefore, the phenology-based GPP estimation substantially improved the estimation in croplands, grasslands, and wetlands compared to original models. Meanwhile, we also noticed that GPP phe-based and GPP VPM estimates lagged behind GPP EC for herbaceous plant in several years, particularly for the grassland site (US-VAR) in 2013. Examining the camera photos of the extracted SOS and EOS revealed no obvious flaws in the extraction process. When tracing back to the results of phenology extraction using the HGAM method at this grassland site, we found that the predefined threshold may have crossed the reconstructed VI three times rather than twice, as in earlier research (Supplementary Materials, Figure S1.19-S1.21, which suggests that the third crossing in the DOY may have been the SOS rather than the first intersection. Then, the third crossing and the second crossing in the DOY were set for the SOS and EOS, and we recalculated the GPP using the same procedures as mentioned above. The results showed that the R 2 between the EVI-based GPP phe-based and GPP EC was not substantially promoted (0.04 vs. 0.02). Therefore, it was still a challenge to use the phenology-based method in GPP estimation for the grassland ecosystem due to its greatly divergent or sharp waves in its VI time series after filtering [18]. Another factor that affects phenology extraction is extreme events (e.g., drought) under Mediterranean climates, resulting in vegetation activity that is significantly higher in winter than in summer, which generates an inflection point which may be potentially misidentified as the SOS or EOS [18]. Additional work should focus on using more grassland sites from additional years to provide a more comprehensive comparison of different phenology methods. The method used to estimate ε0 is another limitation that should not be overlooked. Fitting light response equations is one of the commonly used methods. Specifically, ε0 is estimated from the Michaelis-Menten light response equation to fit observed NEE, PAR, and ecosystem respiration, which was adopted in this study. This method is largely de- The method used to estimate ε 0 is another limitation that should not be overlooked. Fitting light response equations is one of the commonly used methods. Specifically, ε 0 is estimated from the Michaelis-Menten light response equation to fit observed NEE, PAR, and ecosystem respiration, which was adopted in this study. This method is largely dependent on the choice of equation forms. Gan et al. [67] found great variations of ε 0 when comparing its estimates from low light regression and the light response equation across eight biomes. Future work is needed to synthesize these methods to offer a comprehensive vision of ε 0 estimation.

Conclusions
In this study, we investigated and evaluated the performance and robustness of DLF and the HGAM in vegetation phenology extraction from NDVI/EVI/GCC time series across seven PFTs and proposed a phenology-based framework to estimate ε 0 for improved GPP estimation in the LUE model. The results showed that using the HGAM method and EVI data was suited for phenology extraction in forest and cropland ecosystems; wetland phenology extraction preferred the integration of camera photos and the HGAM method; and grassland phenology favored EVI and the DLF method. Based on land surface phenology and light-response curves, we estimated phenology-based ε 0 for the different ecosystems. The results showed that the ranges of ε 0 for different ecosystems were divergent. However, we found high consistency in the seasonal dynamics of ε 0 across these ecosystems. Based on the point, four phenology-based GPP products were estimated from phenology-based ε 0 and four different forms of FAPAR and compared with VPM GPP. In different types of ecosystems, the improved model showed varying degrees of improvement compared to the original model. Generally, the good agreement between the GPP phe-based and GPP EC in our study implies that the maximum LUE estimated in the phenological stages reduces the uncertainty of the parameter to a certain extent and provides a substantial improvement in GPP estimations for cropland, grassland, and wetland.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs15164002/s1, Figure S1.1-S1.21 presents the key phenological stages extracted by HGAM and DLF at the seven flux tower sites. Figure S2.1-S2.21 presents the phenological camera photographs corresponding to the key phenological stages at the seven flux tower sites. Figure S3 presents the linear correlation between the flux tower observation GPP and the four FAPAR proxies.  Data Availability Statement: The flux site observations are from FLUXNET, the phenology camera data are from PhenoCam, and other data presented in this study are available on request from the corresponding author.