Understanding the Land Surface Phenology and Gross Primary Production of Sugarcane Plantations by Eddy Flux Measurements, MODIS Images, and Data-Driven Models

Sugarcane (complex hybrids of Saccharum spp., C4 plant) croplands provide cane stalk feedstock for sugar and biofuel (ethanol) production. It is critical for us to analyze the phenology and gross primary production (GPP) of sugarcane croplands, which would help us to better understand and monitor the sugarcane growing condition and the carbon cycle. In this study, we combined the data from two sugarcane EC flux tower sites in Brazil and the USA, images from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, and data-driven models to study the phenology and GPP of sugarcane croplands. The seasonal dynamics of climate, vegetation indices from MODIS images, and GPP from two sugarcane flux tower sites (GPPEC) reveal the temporal consistency in sugarcane phenology (crop calendar: green-up dates and harvesting dates) as estimated by the vegetation indices and GPPEC data. The Land Surface Water Index (LSWI) is shown to be useful to delineate the phenology of sugarcane croplands. The relationship between the sugarcane GPPEC and the Enhanced Vegetation Index (EVI) is stronger than the relationship between the GPPEC and the Normalized Difference Vegetation Index (NDVI). We ran the Vegetation Photosynthesis Model (VPM), which uses the light use efficiency (LUE) concept and is driven by climate data and MODIS images, to estimate the daily GPP at the two sugarcane sites (GPPVPM). The seasonal dynamics of the GPPVPM and GPPEC at the two sites agreed reasonably well with each other, which indicates that VPM is a powerful tool for estimating the GPP of sugarcane croplands in Brazil and the USA. This study clearly highlights the potential of combining eddy covariance technology, satellite-based remote sensing technology, and data-driven models for better understanding and monitoring the phenology and GPP of sugarcane croplands under different climate and management practices.


Introduction
Sugarcane (complex hybrids of Saccharum spp.) is one of the major cash crops in the world [1]. It accounts for 70% of the world's sugar production [2] and is also a major source of biomass feedstock The sugarcane plantation EC flux tower site (21.6366°S, 47.7897°W) is located at Luiz Antonio municipality, Saõ Paulo State, Brazil [18]. Sugarcane varieties (SP81-3250, SP83-2847, and RB86-7515) were planted in this site over years. The site is 552 m high above sea level and its soil type is sandy soil. The soil and sugarcane stalk yields were representative of southern Brazil [48,49]. The canopy height reached approximately 5 m when sugarcane was at the peak of its growing season. The distance between the sugarcane planting rows was 1.4 m. The sugarcane plantation field (>400 ha) has a gentle slope of <2%, and the close-by vegetation includes pasture, citrus fruit orchards, and savanna woodland [49]. Sugarcane is a multi-year ratoon crop, and on average sugarcane fields in the region have a cropping cycle comprising one plant crop and four ratoon (regrowth) crops [18,49]. At the USR site, the sugarcane plantation was planted in 2003 and replanted in 2007, one year less than the average because a low cane stalk yield was observed in the 2007 harvest (day of the year 149). The detailed agricultural practices can be found in the publication [18].

The Sugarcane Plantation Site in Louisiana, USA
The sugarcane plantation site (29.6341°N, 90.8349°W) is located at the Ardoyne Farm in Schriever, Louisiana, USA, and is managed by the USDA-ARS Sugarcane Research Unit. The HoCP 04-83' variety (reg. no. CV-181, PI 687221) of sugarcane plantation was selected in this site. This site (100 ha) has a long-term (50 yr) history of continuous sugarcane production. The study site is bordered to the north, east, south, and west by 1200, 1320, 360, and 2000 m of continuous sugarcane production fields. The field is graded to 0.2% toward the south. Sugarcane plant rows are spaced at 1.83 m apart. The soil type is Cancienne silty clay loam (Fluvaquentic Epiaquepts). The sugarcane plants usually green up in April and are harvested by December ( Figure S3) [50].  The sugarcane plantation EC flux tower site (21.6366 • S, 47.7897 • W) is located at Luiz Antonio municipality, São Paulo State, Brazil [18]. Sugarcane varieties (SP81-3250, SP83-2847, and RB86-7515) were planted in this site over years. The site is 552 m high above sea level and its soil type is sandy soil. The soil and sugarcane stalk yields were representative of southern Brazil [48,49]. The canopy height reached approximately 5 m when sugarcane was at the peak of its growing season. The distance between the sugarcane planting rows was 1.4 m. The sugarcane plantation field (>400 ha) has a gentle slope of <2%, and the close-by vegetation includes pasture, citrus fruit orchards, and savanna woodland [49]. Sugarcane is a multi-year ratoon crop, and on average sugarcane fields in the region have a cropping cycle comprising one plant crop and four ratoon (regrowth) crops [18,49]. At the USR site, the sugarcane plantation was planted in 2003 and replanted in 2007, one year less than the average because a low cane stalk yield was observed in the 2007 harvest (day of the year 149). The detailed agricultural practices can be found in the publication [18].

The Sugarcane Plantation Site in Louisiana, USA
The sugarcane plantation site (29.6341 • N, 90.8349 • W) is located at the Ardoyne Farm in Schriever, Louisiana, USA, and is managed by the USDA-ARS Sugarcane Research Unit. The HoCP 04-83' variety (reg. no. CV-181, PI 687221) of sugarcane plantation was selected in this site. This site (100 ha) has a long-term (50 yr) history of continuous sugarcane production. The study site is bordered to the north, east, south, and west by 1200, 1320, 360, and 2000 m of continuous sugarcane production fields. The field is graded to 0.2% toward the south. Sugarcane plant rows are spaced at 1.83 m apart. The soil type is Cancienne silty clay loam (Fluvaquentic Epiaquepts). The sugarcane plants usually green up in April and are harvested by December ( Figure S3) [ At the site, a closed-path instrument was used to measure the CO 2 flux. The EC flux tower's height is 9 m. A detailed description of the EC flux tower measurements was given in a previous publication [18]. The GPP was calculated as the difference between the observed NEE and the ecosystem respiration (ER) during daylight hours, which was estimated by a model that uses the exponential relationships between air temperature and nighttime CO 2 fluxes (ER) [51]. We aggregated the half-hourly photosynthetically active radiation (PAR) and CO 2 flux data into the daily sum data. We further averaged the daily climate and carbon flux data to 8-day period data. In this study, the climate and CO 2 flux data of the whole sugarcane growing season in 2005-2007 at this site was used. The daily daytime mean air temperature was defined as the average temperature over the period that has PAR values larger than 10 µmol m −2 s −1 .

The Sugarcane Plantation Site in Louisiana, USA
At this site, an integrated open-path infrared gas analyzer was instrumented with the tower [16]. The calculation of the GPP, nighttime ER, and daytime ER with the fitted exponential equations followed the same procedure as the previous study [51]. We used a similar method to aggregate the daily and 8-day period data from the half-hour data, as was done for the Brazil site. In this study, we used the data of the entire sugarcane growing season in 2017 at this site.
The MOD09A1 Collection 6 product provides 8-day estimates of surface reflectance of seven spectral bands at 500 m spatial resolution. A detailed description of MOD09A1 is given at https: //lpdaac.usgs.gov/products/mod09a1v006/. The MOD09A1 data are available in the Google Earth Engine (GEE), a powerful satellite image editing platform [52,53]. After the data quality control of the MOD09A1 files, we used land surface reflectance data to calculate the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Land Surface Water Index (LSWI) (Equations (1)- (3)). The algorithm documented in the previous publication [54] was used to gap fill the bad quality NDVI and EVI observations. We used times series MODIS data for these two sugarcane sites during 2000-2018, based on the geographical sites (latitude and longitude) of the sugarcane sites. The seasonal and interannual variation of three vegetation indexes during 2000-2018 for these two sugarcane sites is shown in Figure S4.

NCEP Climate Data
The National Centers for Environmental Prediction (NCEP) -Department of Energy (DOE) NCEP/DOE Reanalysis-2 climate dataset was used in our global simulation of VPM model published in Scientific Data paper [55]. The daily climate variables (air temperature, precipitation, and downward shortwave radiation) were downscaled to 500 m using a nonlinear spatial interpolation method [56]. In this study, we used the same climate data for the global VPM simulations and used the daily minimum, mean, and maximum air temperature ( • C), precipitation (mm), and downward shortwave radiation (W m −2 ) for the pixels covering the two sugarcane sites. The 8-day climate data were obtained by calculating the averages of the minimum, mean, and maximum air temperature, precipitation, and downward shortwave radiation over an 8-day period, which is the same temporal resolution as the MODIS data product. A previous study [57] showed a comparison between the NCEP radiation data with in situ radiation measurements at 37 AmeriFlux sites and reported a bias correction factor of 0.8. In this study, we applied this factor to adjust the radiation data. Many studies have indicated that PAR is just a proportion of the SR, and the conversion from SR to PAR ranges from 0.45 to 0.50 [58][59][60].
In this study, we use the mean PAR/SR ratio (0.48) as the conversion ratio [61]. The seasonal dynamics and interannual variation in air temperature ( • C), precipitation (mm day −1 ), and PAR (mol m −2 day −1 ) during 2000-2018was shown in Figure S5.

VPM Model
The VPM model estimates the daily GPP by multiplying the amount of PAR absorbed by chlorophyll in the crop canopy (APAR chl ) and the light use efficiency [37,38]. The model equations are the following: The fraction of PAR absorbed by chlorophyll in the canopy (FPAR chl ) is estimated as a linear function of EVI [55]: The effects of air temperature and water on the maximum light use efficiency (ε 0 ) are calculated as: The parameters ε g and ε 0 refer to the light use efficiency and maximum light use efficiency (µmol CO 2 µmol −1 photosynthetic photon flux density (PPFD) or g C mol −1 PPFD), respectively. The W scalar and T scalar scalars represent the effects of air temperature and water on the light use efficiency of vegetation, respectively. The maximum light use efficiency (LUE) parameter values are obtained from a literature survey and/or analysis of NEE, GPP, and PAR data from CO 2 EC flux tower sites [62]. In this study, we used 0.075 mol of CO 2 mol −1 PPFD (0.9 g C mol −1 PPFD) as the parameter ε 0 for sugarcane (C4 plant), based on the findings of an early study of the CO 2 uptake and quantum yield of photosynthesis in different sugarcane clones at leaf scale [63]. In that study, the authors investigated how leaf nitrogen content affects the quantum yield at the leaf scale for a variety of sugarcane plants, and they found that quantum yield is a linear function of leaf nitrogen content and reached 0.075 mol of CO 2 mol −1 PPFD when the leaf nitrogen content is the largest value in the experiment. In this study, we assume that this quantum yield value at the leaf scale is the maximum LUE parameter at the ecosystem scale.
We calculated the effect of the air temperature on photosynthesis (T scalar ) using Equation (6): In Equation (8), T min , T opt , and T max are the minimum, optimal, and maximum air temperature of photosynthesis, respectively. The T min and T max parameters for these two sugarcane sites were set to be −1 • C and 48 • C, using the same parameter values of the cropland biome in the publication [55]. The site-specific optimum air temperature (T opt ) parameter was estimated by the relationship between the GPP EC or vegetation index (NDVI, EVI) and the daily mean air temperature and daytime mean air temperature. This procedure was already used to estimate the optimum air temperature parameter in studying the carbon fluxes of paddy rice fields [64] and grasslands [65]. According to the results in Section 3.2, the T opt parameter value is approximately~25 • C (both GPP-based and VI-based) at the Brazil sites and~28 • C (GPP-based) or~25 • C (VI-based) at the USA site.
We calculated the effect of water on photosynthesis (W scalar ) and this was estimated by the LSWI, and the maximum LSWI during the plant growing season is assumed to be the LSWI max:

VPM Simulations with the Climate Data from the EC Flux Tower Sites
We used the climate data from the EC flux tower sites to run the VPM simulations (GPP VPM-Site ). The half-hourly air temperature and PAR data were aggregated to daily and 8-day datasets for the model simulation. The resultant GPP VPM-Site data were compared with the GPP EC data.

VPM Simulations with Climate Data from the NCEP Dataset during 2000-2018
We used the NCEP reanalysis 2 dataset to run VPM simulations (GPP VPM-NCEP ) [55]; see Section 2.4 for more details. We calculated and analyzed the resultant GPP data (GPP VPM-NCEP ) and quantified the effect of weather on the GPP of sugarcane croplands during 2000-2018 at these two sugarcane sites.

MODIS GPP and NPP Data Product (MOD17)
The MOD17A2 Gross and Net Primary Productivity (GPP/NPP) product provides daily GPP estimates, and the product uses the light-use efficiency model (PSN) [56]. A description of the MOD17A2 product is given in the website https://lpdaac.usgs.gov/products/mod17a2hv006/. We used the time series GPP data (GPP MOD17A2 ) for these two sugarcane sites during 2000-2018. In this study, we compared the GPP MOD17A2 with the GPP EC , GPP VPM-Site , and GPP VPM-NCEP . The seasonal dynamics and interannual variation of the GPP MOD17A2 during 2000-2018 for two sugarcane sites is shown in Figure S6. Figure 2 shows that the seasonal dynamics of the daily mean air temperature, PAR, and rainfall at the two sugarcane sites over several years. The Brazil site has wet and dry seasons, and the daily mean air temperature is above 15 • C over the year (Figure 2a). The USA site has a temperate climate and a moderately cold winter ( Figure 2b). The seasonal dynamics of the PAR, air temperature, and rainfall were similar, with the peaks in mid-summer at the site ( Figure 2b).

The Seasonal Dyanmics of Climate, Vegetation Indices, and GPP at the Two Sugarcane Tower Sites
The land surface phenology (LSP) includes the start of the growing season (SOS), the end of growing season (EOS), and the growing season length (GSL). The seasonal dynamics of the vegetation indices (NDVI, EVI, and LSWI) at the two sites reveals the LSP metrics in terms of the canopy structure and the process of crop development during the growth season (Figure 3a (Figure 3a). At the USA site, the vegetation indices started to rise rapidly in late April, reached their maximum value by August, and gradually dropped to <0.4, <0.1, and <0 by December 2017, respectively ( Figure 3b). Using the same criteria (EVI > 0.1, LSWI > 0; EVI < 0.1, LSWI < 0), the VI-based SOS and EOS of sugarcane at the USA site were April and December, respectively, and the GSL was about 8 months. The differences in the VI-based SOS VI , EOS VI , and GSL VI between these two sites are largely determined by the air temperature and management practice (e.g., harvest time) (Figure 2).  (Figure 3a). At the USA site, the vegetation indices started to rise rapidly in late April, reached their maximum value by August, and gradually dropped to < 0.4, < 0.1, and < 0 by December 2017, respectively ( Figure 3b). Using the same criteria (EVI > 0.1, LSWI > 0; EVI < 0.1, LSWI < 0), the VI-based SOS and EOS of sugarcane at the USA site were April and December, respectively, and the GSL was about 8 months. The differences in the VI-based SOSVI, EOSVI, and GSLVI between these two sites are largely determined by the air temperature and management practice (e.g., harvest time) ( Figure 2).
The seasonal dynamics of NEEEC and GPPEC at the two sites reveal the phenology of sugarcane croplands from the perspective of ecosystem function (Figure 3c, d). At the Brazil site, the GPPEC rose steadily in May 2005, when the GPPEC was larger than 1 g C m -2 day -1 ; it had a maximum level in December 2005 and decreased to zero by May 2006 (Figure 3c). Using the criterion of GPP >= 1 g C/m 2 /day, the GPP-based SOS of the sugarcane was in mid-May in 2005 and the EOS was in May 2006, which corresponds well with the growing season delineated by the vegetation indices. The GPPEC started to increase again after the onset of the rainy season in October 2006 [18] and reached the maximum in February of 2007 (Figure 3c). At the USA site, the GPPEC rose rapidly in April (GPPEC larger than 1 g C m −2 day −1 ), had a maximum value in August, and decreased to zero by December (Figure 3d). Using the criterion of GPP ≥ 1 g C/m 2 /day, the GPP-based sugarcane growing season at  (Figure 3c). At the USA site, the GPP EC rose rapidly in April (GPP EC larger than 1 g C m −2 day −1 ), had a maximum value in August, and decreased to zero by December (Figure 3d). Using the criterion of GPP ≥ 1 g C/m 2 /day, the GPP-based sugarcane growing season at the USA site ranged from April (SOS) to December (EOS) (Figure 3d), which corresponds well with the growing season by the vegetation indices (Figure 3b). The growing season length (GSL) of the sugarcane, defined as the number of days with a GPP EC larger than 1 g C m −2 day −1 , lasted approximately 12 months at the Brazil site and 8months at the USA site. The differences in the GPP-based SOS GPP , EOS GPP , and GSL GPP between these two sites are largely determined by the air temperature and management practice (e.g., harvest time) (Figure 2)

Estimation of the Site-Specific Optimum Air Tmeperature (T opt ) for GPP at the Two Sugarcane Sites
We assessed the biophysical performance of the vegetation indices at the sugarcane sites in terms of the GPP dynamics. Figure 4 shows the correlation between the GPP EC and vegetation indices (NDVI and EVI) during the sugarcane growing seasons at the two sites. The GPP EC has slightly stronger linear relationships with EVI (R 2 = 0.68 at the Brazil site and R 2 = 0.79 at the USA site) than with the NDVI (R 2 = 0.52 at the Brazil site and R 2 = 0.76 at the USA site) (Figure 4). The slightly lower R 2 values at the Brazil sugarcane sites is in part attributed to the effect of frequent cloud cover and shadow on vegetation indices during the rainy season, as evidenced by the relatively large variation in the vegetation indices (Figure 2a).
Remote Sens. 2020, 12, 8 of 22 the USA site ranged from April (SOS) to December (EOS) (Figure 3d), which corresponds well with the growing season by the vegetation indices (Figure 3b). The growing season length (GSL) of the sugarcane, defined as the number of days with a GPPEC larger than 1 g C m −2 day −1 , lasted approximately 12 months at the Brazil site and 8months at the USA site. The differences in the GPPbased SOSGPP, EOSGPP, and GSLGPP between these two sites are largely determined by the air temperature and management practice (e.g., harvest time) (Figure 2)

Estimation of the Site-Specific Optimum Air Tmeperature (Topt) for GPP at the Two Sugarcane Sites
We assessed the biophysical performance of the vegetation indices at the sugarcane sites in terms of the GPP dynamics. Figure 4 shows the correlation between the GPPEC and vegetation indices (NDVI and EVI) during the sugarcane growing seasons at the two sites. The GPPEC has slightly stronger linear relationships with EVI (R 2 = 0.68 at the Brazil site and R 2 = 0.79 at the USA site) than with the NDVI (R 2 =0.52 at the Brazil site and R 2 =0.76 at the USA site) (Figure 4). The slightly lower R 2 values at the Brazil sugarcane sites is in part attributed to the effect of frequent cloud cover and shadow on vegetation indices during the rainy season, as evidenced by the relatively large variation in the vegetation indices (Figure 2a).  We evaluated the effects of air temperature on the GPPEC by the relationships between the GPPEC and the daily daytime mean air temperature and daily mean air temperature at the two sites ( Figure  5). At the Brazil site, the GPPEC increased as the daily daytime mean air temperature or daily mean air temperature increased and reached its plateau at ~25 °C (Figure 5a,c). At the USA site, the GPPEC increased as the daily daytime mean air temperature or daily mean air temperature increased and reached its plateau at 28 °C (Figure 5b,d). We evaluated the effects of air temperature on the GPP EC by the relationships between the GPP EC and the daily daytime mean air temperature and daily mean air temperature at the two sites ( Figure 5). At the Brazil site, the GPP EC increased as the daily daytime mean air temperature or daily mean air temperature increased and reached its plateau at~25 • C (Figure 5a,c). At the USA site, the GPP EC increased as the daily daytime mean air temperature or daily mean air temperature increased and reached its plateau at 28 • C (Figure 5b,d). We further investigated the effect of air temperature (daily daytime mean air temperature, daily mean air temperature) on the vegetation indices (NDVI and EVI) ( Figure 6). The NDVI and EVI rose with the daily daytime mean air temperature and reached their plateau at 25 °C at the Brazil location and 25 °C for the USA location (Figure 6a, b). The NDVI and EVI were positively correlated with the daily mean air temperature and reached their plateau at 25 °C at the Brazil site and 28°C at the USA site (Figure 6c, d).
For simulations of the VPM model (see Section 3.3), the optimum air temperature (Topt) at the Brazil site is set at 25 °C (both GPP-based and EVI-based) and the optimum air temperature (Topt) at the USA site is set 28 ° C (GPP-based) or 25 °C (EVI-based), respectively, as part of a model sensitivity test. We further investigated the effect of air temperature (daily daytime mean air temperature, daily mean air temperature) on the vegetation indices (NDVI and EVI) ( Figure 6). The NDVI and EVI rose with the daily daytime mean air temperature and reached their plateau at 25 • C at the Brazil location and 25 • C for the USA location (Figure 6a,b). The NDVI and EVI were positively correlated with the daily mean air temperature and reached their plateau at 25 • C at the Brazil site and 28 • C at the USA site (Figure 6c,d).
For simulations of the VPM model (see Section 3.3), the optimum air temperature (T opt ) at the Brazil site is set at 25 • C (both GPP-based and EVI-based) and the optimum air temperature (T opt ) at the USA site is set 28 • C (GPP-based) or 25 • C (EVI-based), respectively, as part of a model sensitivity test.

Seasonal Dynamics of GPP as Simulated by the VPM Model with the Climate Data from the EC Flux Tower Sites (GPP VPM-Site )
The seasonal dynamics of the GPP EC and GPP VPM-Site at the two sites agreed reasonably well with each other (Figure 7 (Figure 7a). At the USA site, the GPP VPM-Site increased from April, reached maximum its value at August, and dropped by late December. The maximum daily G PPVPM-Site was slightly lower than the GPP EC for July-August (Figure 7b).
There are strong linear relationships between the GPP EC and GPP VPM-Site in the sugarcane growing seasons at the two sites (Figure 7c and 1976 g C m −2 yr −1 (28 • C), or 10.06~10.17% lower than the sum of the GPP EC (2200 g C m −2 yr −1 ). The seasonal sum of the GPP VPM-Site during the growing season defined by LSWI at the USA site in 2017 was between 1699 g C m −2 yr −1 (25 • C) and 1709 g C m −2 yr −1 (28 • C), or 9.68~10.25% lower than the sum of the GPP EC (1893 g C m −2 yr −1 ) ( Table 1).

Seasonal Dynamics of GPP as Simulated by the VPM Model with the Climate Data from the EC Flux Tower Sites (GPPVPM-Site)
The seasonal dynamics of the GPPEC and GPPVPM-Site at the two sites agreed reasonably well with each other (Figure 7 (Figure 7a). At the USA site, the GPPVPM-Site increased from April, reached maximum its value at August, and dropped by late December. The maximum daily GPPVPM-Site was slightly lower than the GPPEC for July-August (Figure 7b).
There are strong linear relationships between the GPPEC and GPPVPM-Site in the sugarcane growing seasons at the two sites (Figure 7c (Table 1).

Interannual Variation in GPP as Simulated by the VPM Model with NCEP Climate Data (GPPVPM-NCEP)
The seasonal dynamics and interannual variation in GPPVPM-NCEP from the VPM model with the NCEP climate data at the two sites during 2000-2018 reflects the effect of climate (drought) on the sugarcane GPP (Figure 8). We calculate the annual GPPVPM-NECP of sugarcane plantations as the sum of all the GPP data with GPPVPM≥1 g Cm

Biophysical Performance of Vegetation Indices for the Sugarcane Plantations
Satellite-based vegetation indices (VIs), including greenness-associated VIs (e.g., NDVI, EVI) and water-associated VIs (e.g., NDWI or LSWI) [66][67][68], have been used as proxies for several biochemical and biophysical variables including leaf area index, canopy chlorophyll content, and gross primary production [69]. Our results show that at the sugarcane croplands, the linear relationship between the GPPEC and EVI is slightly stronger than that between the GPPEC and NDVI. Similar results were reported in other studies of paddy rice fields [44], maize croplands [42], grasslands [41], and forests [36][37][38].
Crop phenology or calendar information (e.g., planting date, harvest date) are useful for crop management, crop yield estimation, and carbon cycle study [70,71]. Several studies used reflectance and vegetation indices data to track or classify the sugarcane area over several years and predict the cane stalk yield [72][73][74][75][76][77]. Time series MODIS-derived EVI data were used to classify sugarcane croplands in Brazil [72]. The land surface phenology of croplands can be tracked and delineated by the seasonal dynamics of vegetation indices (NDVI, EVI, and/or LSWI) from the perspective of ecosystem structure (leaf area index, greenness) [38,39]. Our study shows that NDVI, EVI, and LSWI track the phenological dynamics of sugarcane croplands. According to the detailed information on sugarcane cultivation history at the two sites [16,18,50,78], the fields have been used to grow sugarcane for at least a few decades (Table S1). The time series vegetation index data from MODIS during 2000-2018 ( Figure S4) show that the LSWI dropped to < 0 when the sugarcane fields were at the harvest stage at both the sites, which clearly showed that in the fall/winter season the date after the first negative LSWI values (< 0) corresponds to the harvesting date or the ending date in the plant growing season. As shown in Figure 3 and Figure S4, the LSWI-based algorithms accurately identified the harvest dates at the Brazil site and the USA site. This is consistent with our previous studies that used the LSWI to delineate phenological metrics [42,44,79]. In addition, the land surface

Biophysical Performance of Vegetation Indices for the Sugarcane Plantations
Satellite-based vegetation indices (VIs), including greenness-associated VIs (e.g., NDVI, EVI) and water-associated VIs (e.g., NDWI or LSWI) [66][67][68], have been used as proxies for several biochemical and biophysical variables including leaf area index, canopy chlorophyll content, and gross primary production [69]. Our results show that at the sugarcane croplands, the linear relationship between the GPP EC and EVI is slightly stronger than that between the GPP EC and NDVI. Similar results were reported in other studies of paddy rice fields [44], maize croplands [42], grasslands [41], and forests [36][37][38].
Crop phenology or calendar information (e.g., planting date, harvest date) are useful for crop management, crop yield estimation, and carbon cycle study [70,71]. Several studies used reflectance and vegetation indices data to track or classify the sugarcane area over several years and predict the cane stalk yield [72][73][74][75][76][77]. Time series MODIS-derived EVI data were used to classify sugarcane croplands in Brazil [72]. The land surface phenology of croplands can be tracked and delineated by the seasonal dynamics of vegetation indices (NDVI, EVI, and/or LSWI) from the perspective of ecosystem structure (leaf area index, greenness) [38,39]. Our study shows that NDVI, EVI, and LSWI track the phenological dynamics of sugarcane croplands. According to the detailed information on sugarcane cultivation history at the two sites [16,18,50,78], the fields have been used to grow sugarcane for at least a few decades (Table S1). The time series vegetation index data from MODIS during 2000-2018 ( Figure S4) show that the LSWI dropped to <0 when the sugarcane fields were at the harvest stage at both the sites, which clearly showed that in the fall/winter season the date after the first negative LSWI values (<0) corresponds to the harvesting date or the ending date in the plant growing season. As shown in Figure 3 and Figure S4, the LSWI-based algorithms accurately identified the harvest dates at the Brazil site and the USA site. This is consistent with our previous studies that used the LSWI to delineate phenological metrics [42,44,79]. In addition, the land surface phenology of croplands can also be tracked and delineated by the seasonal dynamic of GPP EC and NEE EC data at the EC flux tower sites from the perspective of ecosystem function (physiology) [36,37,41]. Our study shows that GPP EC tracks well the carbon uptake period using a simple criterion (GPP EC ≥ 1 g C /m 2 /day). It is important to note that this study demonstrates the temporal consistency in the land surface phenology derived from these two approaches. The simulations of the LUE-based GPP models driven by in situ climate and known C3 or C4 maximum LUE parameters were assessed at multiple crop EC flux tower sites [44,45,[80][81][82]. The results of this study show that our VPM model has a good performance in modelling the GPP of sugarcane croplands when using a C4 maximum LUE parameter, site-specific T opt parameter, and local climate data (GPP VPM-Site ) or global climate data (GPP VPM-NCEP ). Several global GPP datasets derived from satellite-based LUE models are also available to the community, including the Terra/MODIS Gross and Net Primary Production (GPP/NPP) data product (MOD17A2) [56,83] and the GPP dataset from VPM simulations (GPP VPM-globe ) [55]. We compared the seasonal dynamic of the GPP EC , GPP VPM-Site , GPP VPM-NCEP , GPP MOD17A2 , and GPP VPM-Globe at the two sites ( Figure 10). The GPP MOD17A2 was much lower than the GPP EC , GPP VPM-Site , and GPP VPM-Globe during the plant growing season at the two sites (Figure 10a,b). The scatterplots of GPP MOD17A2 , GPP VPM-Site , GPP VPM-NCEP , and GPP VPM-Globe versus the GPP EC over the plant growing seasons at the two sites showed strong linear correlations between the GPP EC and other GPP estimates (Figure 10c,d). Using the GPP EC as reference, the GPP MOD17A2 underestimated the GPP by 61% at the USA site and 59% at the Brazil site. The GPP VPM-Globe underestimated the GPP by 35% at the USA site and 40% at the Brazil site. The large underestimation by the global data products (GPP MOD17A2 and GPP VPM-Globe ) at the two sites can be attributed to the maximum LUE parameter used in the global simulations. For the MOD17A2 GPP/NPP product, the C3 maximum LUE parameter was used in its model simulation, as the land cover map using (MCD12C1 C55) does not separate the C3 and C4 crops. For the GPP VPM-Globe product [55], the VPM simulations used a mix of C3 and C4 maximum LUE parameters for croplands and the land cover maps (MCD12C1, C55). The maximum LUE parameter values of the C4 plants are substantially higher than those of the C3 plants [84]. The results from this study clearly highlight the urgent need for developing annual maps of C3 and C4 crops so that model simulations can use the appropriate C3 or C4 maximum LUE parameter (ε0) for specific crop function types. Such comparison analyses also suggest that the global GPP data products to estimate the sugarcane GPP at local and regional scales should be used cautiously.

Sources of Errors and Uncertainty in Predicted GPP VPM-site at Sugarcane EC Flux Tower Sites
VPM simulations at the two sugarcane sites have several sources of error or uncertainty, including in situ climate datasets, time series MODIS vegetation indices, model parameters, and land cover types within the MODIS image pixels and the footprint of EC flux tower sites. The maximum LUE parameter values are estimated by various estimation methods and directly affect GPP simulations at the ecosystem level [85,86]. Several previous in situ studies have reported that the LUE varied under rain-fed and irrigated conditions-for example, for sugarcane varieties in Sri Lanka, the LUE parameter ranged between 1.63 and 2.09 g dry biomass MJ -1 under irrigated and between 0.71 and 1.03 g dry biomass MJ −1 under rain-fed conditions [87]. Another publication reported that the maximum LUE for sugarcane at tropical region in Australia ranges from 1.7 g dry biomass MJ −1 [88] to 2.0 g dry biomass MJ −1 [89]. One field experiment in Brazil showed that maximum LUE parameter ranged between 1.74 and 2.28 g C MJ −1 (0.85-1.12 g C mol −1 PPFD) [90]. In the above-mentioned papers, the LUE was calculated using the linear regression between the cumulative radiation intercepted and the aboveground dry weight. In addition, many studies analyzed the PAR and gross primary production data to estimate the maximum LUE parameter (LUE GPP ) [21]. In the above publication, the LUE AGB or LUE GPP parameter at two sugarcane sites in Hawaii, USA, were 0.88 ± 0.02 g C MJ −1 (0.43 ± 0.01 g C mol −1 PPFD) or 1.24 ± 0.22 g C MJ −1 (0.61 ± 0.11 g C mol −1 PPFD) at the site located at a low altitude and 0.75 ± 0.01 g C MJ −1 (0.37 ± 0.01 g Cmol −1 PPFD) or 1.15 ± 0.15 g C MJ −1 (0.56 ± 0.07 g Cmol −1 PPFD) at the site located at a high altitude. Different LUE definitions should be considered a source of errors associated with satellite-based GPP estimates. In our study, the maximum LUE parameter at the Brazil site and the USA site was set to 0.9 g Cmol −1 PPFD, which was slightly lower than the maize crop (0.92 g Cmol −1 PPFD), another C4 crop in China [79]. The range value of the maximum LUE parameter for sugarcane cropland indicates the need for the evaluation of the maximum LUE parameter over diverse sugarcane croplands in the world under different climate conditions and management practices.

Sources of Errors and Uncertainty in Predicted GPPVPM-site at Sugarcane EC Flux Tower Sites
VPM simulations at the two sugarcane sites have several sources of error or uncertainty, including in situ climate datasets, time series MODIS vegetation indices, model parameters, and land cover types within the MODIS image pixels and the footprint of EC flux tower sites. The maximum LUE parameter values are estimated by various estimation methods and directly affect GPP simulations at the ecosystem level [85,86]. Several previous in situ studies have reported that the LUE varied under rain-fed and irrigated conditions-for example, for sugarcane varieties in Sri Lanka, the LUE parameter ranged between 1.63 and 2.09 g dry biomass MJ -1 under irrigated and between 0.71 and 1.03 g dry biomass MJ −1 under rain-fed conditions [87]. Another publication reported that the maximum LUE for sugarcane at tropical region in Australia ranges from 1.7 g dry biomass MJ −1 [88] to 2.0 g dry biomass MJ −1 [89]. One field experiment in Brazil showed that maximum LUE parameter ranged between 1.74 and 2.28 g C MJ −1 (0.85-1.12 g C mol −1 PPFD) [90]. In the above-mentioned papers, the LUE was calculated using the linear regression between the cumulative radiation intercepted and the aboveground dry weight. In addition, many studies analyzed the PAR and gross primary production data to estimate the maximum LUE parameter (LUEGPP) [21]. In the above publication, the LUEAGB or LUEGPP parameter at two sugarcane sites in Hawaii, USA, were 0.88 ± 0.02 g C MJ −1 (0.43 ± 0.01 g C mol -1 PPFD) or 1.24 ± 0.22 g C MJ −1 (0.61 ± 0.11 g C mol −1 PPFD) at the site located at a low altitude and 0.75 ± 0.01 g C MJ −1 (0.37 ± 0.01 g Cmol −1 PPFD) or 1.15 ± 0.15 g C MJ −1 (0.56 ± 0.07 g Cmol −1 PPFD) at the site located at a high altitude. Different LUE definitions should be considered a source of errors associated with satellite-based GPP estimates. In our study, the maximum LUE parameter at the Brazil site and the USA site was set to 0.9 g Cmol −1 PPFD, which was

Conclusions
This is the first case study that combines in situ data from sugarcane EC flux towers and satellite (MODIS) images to estimate the GPP of sugarcane croplands using the VPM model. The sensitivity analyses between the GPP EC and vegetation indices (NDVI and EVI) at the two sites provide new insight into the biophysical performance of vegetation indices. The sensitivity analysis between the GPP EC and temperature variables, including the mean daily air temperature, mean daytime air temperature, land surface temperature, and the analysis between the EVI and temperature variables at the two sites provides a reasonable method for finding the optimum air temperature values at individual sites, which may improve VPM model simulations at the landscape scale. The results demonstrate that the VPM model performs well in estimating the GPP seasonal dynamics at the two sugarcane sites. The maximum LUE parameter value in this study (0.9 g C mol −1 PPFD) is reasonable for applying the VPM model at these sugarcane croplands. The linked analyses of the flux data (GPP EC ) from the EC flux tower sites, predicted GPP (GPP VPM-Site ) with in-situ climate data from the EC flux tower site, and predicted GPP (GPP VPM-NCEP ) with the NCEP climate data clearly demonstrates the feasibility of the VPM model to estimate and predict sugarcane physiological parameters in different climate systems. Additional evaluations of the VPM model at other sugarcane EC flux tower sites are still needed, which would contribute to a better understanding of the potential sources of errors and uncertainty in the simulations of the VPM model.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/14/2186/s1, Figure S1: Averaged sugarcane production from the top 10 producer countries in the world (Data source: FAO, averaged values over the period of 1994-2017); Figure S2: Annual harvested area, production and yield of sugarcane in Brazil and USA during1961-2017; Figure S3: Sugarcane production cycle at the Louisiana site, USA; Figure Table S1: A detailed information on sugarcane planting history at the two sugarcane sites during 2000-2018; Table S2: A summary of annual precipitation, mean annual air temperature, annual PAR from NCEP Reanalysis-2 climate data, and GPP estimated from the VPM model with NCEP Reanalyssi-2 climate data (GPP VPM-NCEP ) during the sugarcane plant growing season at the sugarcane flux tower sites in Brazil and USA.