Improving MODIS Gross Primary Productivity by Bridging Big‐Leaf and Two‐Leaf Light Use Efficiency Models

Gross primary productivity (GPP) is an important component of the terrestrial carbon cycle in climate change research. The global GPP product derived using Moderate Resolution Imaging Spectroradiometer (MODIS) data is perhaps the most widely used. Unfortunately, many studies have indicated evident error patterns in the MODIS GPP product. One of the main reasons for this is that the applied big‐leaf (BL) MOD17 model cannot properly handle the variable relative contribution of sunlit and shaded leaves to the total canopy‐level GPP. In this study, we developed a model for correcting the errors in the MODIS GPP product by bridging BL and two‐leaf (TL) light use efficiency (LUE) models (CTL‐MOD17). With the available MODIS GPP product, which considers environmental stress factors, the CTL‐MOD17 model only needs to reuse the two inputs of the leaf area index (LAI) and incoming radiation. The CTL‐MOD17 model was calibrated and validated at 153 global FLUXNET eddy covariance (EC) sites. The results indicate that the modeled GPP obtained with the correction model matches better with the EC GPP than the original MODIS GPP product at different time scales, with an improvement of 0.07 in R2 and a reduction in root‐mean‐square error (RMSE) of 117.08 g C m−2 year−1. The improvements are more significant in the green season when the contribution of shaded leaves is larger. In terms of the global spatial pattern, the obvious underestimation in the regions with high LAI and the overestimation in the low LAI regions of the MODIS GPP product is effectively corrected by the CTL‐MOD17 model. This paper not only bridges the BL and TL LUE models, but also provides a new and simple method to obtain accurate GPP through reusing two inputs used in producing the MODIS GPP product.


Introduction
Gross primary productivity (GPP) is the total amount of organic carbon fixed by green plants through photosynthesis per unit of space and time.GPP is a principal indicator of biosphere carbon fluxes and plays an essential role in regulating the global terrestrial carbon cycle (Beer et al., 2010;Chen et al., 2019;Marcott et al., 2014;Wei et al., 2023).Therefore, high-accuracy GPP products at regional or global scales are crucial for understanding the changes of terrestrial ecosystem carbon cycling and its response to global change (Kolby Smith et al., 2016;Running et al., 2004;Wang et al., 2020).
In the past decades, numerous methods have been developed to obtain GPP, and a number of relevant products with various accuracies have been published, which mainly include in situ measurement based and model-based products (Hu et al., 2022;Marshall et al., 2018;Tang & Zhuang, 2009;Xie et al., 2020).The eddy covariance (EC)-based GPP products are the most representative in situ measurement based products derived from the measured net ecosystem exchange of carbon dioxide (CO 2 ) (Reichstein et al., 2005;Wutzler et al., 2018).Typical examples are FLUXNET2015, AmeriFlux, EuroFlux, and ChinaFlux (Kumar et al., 2016;Pastorello et al., 2020).These products provide accurate and continuous site-level data, but they are limited by the sparse and small spatial footprints (Barcza et al., 2020).Model-based products, including process-based, light use efficiency (LUE)based, solar-induced chlorophyll fluorescence (SIF)-based, and machine learning (ML)-based products, can describe the regional and global distributions of GPP.The SIF-based and ML-based products are produced using emerging methods and have been widely used in recent years.However, the SIF-based products are limited by the short time range and the spatial discontinuity of satellite observations (Frankenberg et al., 2014;Köhler et al., 2018), and the ML methods cannot explain the mechanism of the photosynthetic process and depend on the data used for the model training (Jung et al., 2011(Jung et al., , 2019;;Zhang et al., 2021).The process-based products are based on models with solid biophysical mechanisms, but the massive inputs and parameter requirements usually result in an insufficient spatial resolution and poor accuracy for these products, such as the CLM5 and BIOME-BGC products (Huang et al., 2011;Lawrence et al., 2019).The LUE model based products include the Moderate Resolution Imaging Spectroradiometer (MODIS) MOD17 product (Running et al., 2004), the Global Land Surface Satellite (GLASS) product (Yuan et al., 2010), and the TL-LUE product (Bi & Zhou, 2022).
The MOD17 product is probably the most widely used GPP product at present, but many studies have indicated significant error patterns in the MOD17 GPP product (He et al., 2013;Zhang et al., 2012).Significant underestimation can be found in areas with high LAI and during the growing season, while overestimation is found in low LAI areas, which greatly impacts the reliability of the subsequent applications (Chen et al., 2020;Ma et al., 2014;Zhang et al., 2012).The sources of the uncertainties come from many aspects, including the unreliable input data, the scale differences of the multi-source data, the complications of the background reflectance, the limitations of the environmental stressors, and the vegetation type-specific parameters of the used model, and so on (Baldocchi & Penuelas, 2019;Pei et al., 2022;Serbin et al., 2013;Wagle et al., 2016).Among the different uncertainties, the inaccuracy of the input data and parameters can have a great impact.For example, many studies have indicated that the MODIS LAI induces large uncertainties in the MOD17 GPP product (Heiskanen et al., 2012;Knyazikhin et al., 1998;Serbin et al., 2013;Yan et al., 2016).However, these issues are also inevitable in other models, and the model structure of the MOD17 model may be the unique source of error in this product (Heiskanen et al., 2012;Knyazikhin et al., 1998;Yan et al., 2021;Yuan et al., 2014).The MOD17 model is a "big-leaf" (BL) LUE model that treats the entire vegetation canopy as a BL and only uses a constant maximum LUE and absorbed photosynthetically active radiation (APAR) value to represent the average canopy conditions of vegetation.Sunlit leaves absorb both direct and diffuse radiation simultaneously, whereas shaded leaves only absorb diffuse radiation, and hence the GPP values of sunlit and shaded leaves differ considerably (Chen et al., 1999;Zhang et al., 2012).The contribution of shaded leaves to the canopy-level GPP can differ greatly, depending on the amount of the LAI.In dense vegetation with high LAI, the shaded leaf area can greatly exceed the sunlit leaf area, resulting in the shaded leaves making a greater contribution to canopy-level GPP than the sunlit leaves, while in sparse vegetation with low LAI, most of the leaves are sunlit, and the sunlit leaves contribute the most to the canopy-level GPP.However, a BL model assumes a canopy to be a BL, without the shaded leaf component, and is inherently incapable of considering the different contributions of shaded and sunlit leaves.A BL model can consider the average level of the shaded leaf contribution by finding an appropriate value of the maximum LUE, but it will underestimate canopy-level GPP when the shaded leaf contribution is more than the average (i.e., high LAI) and will overestimate the canopy-level GPP when the shaded leaf contribution is less than the average (i.e., low LAI) (Zhang et al., 2012;Zhou et al., 2016).
This problem can mostly be solved by the TL LUE models, which are developed by separately calculating the GPP of sunlit and shaded leaves, based on the approach in the Biosphere-atmosphere Exchange Process Simulator (BEPS), which was renamed from the Boreal Ecosystem Productivity Simulator (Alexandre et al., 1997;Chen et al., 1999).Typical examples of the TL LUE models are the TL-LUE model (He et al., 2013), the mountainous two-leaf (MTL-LUE) model (Xie & Li, 2020), the TL LUE model considering a radiation scalar (RTL-LUE) (Guan et al., 2021(Guan et al., , 2022)), and the revised EC-LUE model (Zheng et al., 2019).These models have been shown to outperform BL models across a wide range of biomes by correcting the underestimation within high GPP ranges by considering the contribution of shaded leaves.Although the various TL methods have already been widely used and achieved satisfactory results, they also require different meteorological inputs, which are usually obtained from reanalysis data with a coarse resolution or through interpolation from limited meteorological stations (Xie & Li, 2020;Zhou et al., 2016).Meanwhile, the existing GPP products already include information on environmental stressors, and although some uncertainties also exist, some products even use exactly the same environmental stressors as the related TL LUE models.If we improve the BL MOD17 GPP by only correcting the bias in the canopy APAR based on TL theory, that is, using exactly the same environmental stress functions with the MOD17 model as the previous TL LUE models do, it would be unnecessary to use temperature and moisture related meteorological data to calculate the environmental stressors because this information is already included in the MOD17 product, and the environmental input data requirements can be greatly reduced.
Since there are already many GPP products published, taking these products as a basis and correcting their apparent biases would be an effective and useful approach for producing new high-accuracy GPP products (Liang et al., 2021;Zhang & Ye, 2021;Zhang et al., 2017).The MOD17 GPP data are popularly used as the MODIS official product with well-recognized biases, so it is helpful to improve this product based on a simple correction procedure (Hwang et al., 2008;Li et al., 2021).Most of the previous TL LUE models were developed on the basis of the MOD17 model, so it may be naturally practicable to correct the uncertainties regarding the calculation of canopy APAR in the MOD17 GPP product by integrating the contributions of sunlit and shaded leaves, so as to bridge the gap between the BL and TL LUE models (Guan et al., 2022;He et al., 2013).In this regard, a similar GPP estimation accuracy to the TL models can be achieved based on the existing MODIS product using only a few additional data inputs.Furthermore, although many BL and TL LUE models have already been developed, they are all independent of each other, and it is still unknown whether a conversion model from BL to TL simulation is possible or not.The superiority of TL models over BL models demonstrated in previous studies was mostly established based on separate model simulations.If we can prove that simply correcting a BL product based on TL theory can achieve a better GPP product, the previous conclusions regarding the superiority of TL models could be consolidated (Xie & Li, 2020;Zhou & Yan, 2016).Therefore, it is of great significance to explore the possibility of bridging the gap between these two types of models.
The objectives of this study were: (a) to develop a correction model (CTL-MOD17) to obtain GPP with high accuracy based on the MODIS GPP product; (b) to evaluate the accuracy of the CTL-MOD17 model and analyze the advantage of TL theory at a global scale from a correction perspective; and (c) to investigate the uncertainties in the MOD17 product and the contribution of the used BL model.

FLUXNET Data
EC data from the FLUXNET2015 data set (Baldocchi et al., 2018;Falge et al., 2006;Pastorello et al., 2020) collected from multiple flux networks across the globe (www.fluxdata.org)were used to validate the CTL-MOD17 model.As in many previous studies, GPP partitioned by the daytime (DT) method using a hyperbolic light response curve (termed "GPP_DT_VU-T_REF") was employed as the reference GPP to evaluate the model performance (Papale et al., 2006;Reichstein et al., 2005).In total, 153 sites and 1,191 site-years covering 11 vegetation types (Table S1 and Figure S1 in Supporting Information S1) were selected in this study.The main principles of selection were as follows: (a) more than two years of observations over 2001-2015; and (b) more than 80% valid data in the vegetation green season (April-October in the Northern Hemisphere and the remaining months in the Southern Hemisphere) (Liu et al., 2013;Piao et al., 2007).For each site, tower-based GPP data over each 8-day interval were summed to match the time step of the MODIS GPP product, and the FLUXNET2015 data set was randomly divided into two parts, with 20% of the site-years used for the validation and the remaining data used for the parameter optimization.The number of sites and site-years for each vegetation type are listed in Table 1.

MODIS Data
The MOD17A2H GPP and MOD15A2H LAI products from 2001 to 2015 were used to drive the CTL-MOD17 model (https://ladsweb.modaps.eosdis.nasa.gov/).Both products are 8-day composites with a spatial resolution of tower were extracted to match the footprint of the EC data (Barcza et al., 2020;Zhou et al., 2016).Furthermore, to reduce the impacts of unrealistic short-term fluctuations caused by cloud contamination and other atmospheric impacts, the MODIS LAI time-series data were smoothed using the temporal filtering method (Chen et al., 2004;Chu et al., 2021).

Other Data
In order to agree with the data in the MODIS GPP production scheme, the solar incident shortwave radiation data from the second Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) tavg1_2d_rad_Nx reanalysis data set were also used (https://gmao.gsfc.nasa.gov/).Although uncertainties are indicated by many studies in the average bias of the regionally weighted radiation, it is also used in this study to keep consistent with MOD17 product (Running et al., 2004;Wang et al., 2023;Zhao et al., 2006).The MERRA-2 data set is hourly time-averaged with a spatial resolution of 0.5°× 0.625°.The radiation data in Greenwich Mean Time in this data set were converted to local time at each flux site.The hourly data were then composited into 8day data, and the spatial resolution was interpolated to 500 m to match the spatial resolution of the LAI and MOD17 GPP data.
The TL-LUE GPP data set was also used in this study, in order to compare the global GPP spatial distribution with the MOD17 product and CTL-MOD17 correction results (Bi & Zhou, 2022).The TL-LUE GPP product is produced using an updated TL-LUE model, considering the atmospheric CO 2 concentration regulation scalar and air temperature regulation scalar (Bi & Zhou, 2022).The TL-LUE GPP product has demonstrated consistently high-precision estimates across ecosystems of various vegetation types, especially in areas with a high vegetation density.This product is provided at a 0.05°spatial resolution every 8 days from 1992 to 2020 (https://doi.org/10.5061/dryad.dfn2z352k)(Bi & Zhou, 2022).

Description of the MOD17 Model
The MOD17 algorithm was developed based on Monteith's radiation conversion efficiency concept from 1972.GPP is calculated as (Running et al., 2004): where ε max is the maximum LUE with the different vegetation types; f (VPD) and g (T min ) are the scalars of vapor pressure deficit (VPD) and minimum air temperature, respectively; and f PAR is the fraction of PAR absorbed by the canopy.These are calculated as follows: where VPD max , VPD min , TMIN min , and TMIN max are the meteorological parameters dependent on the vegetation type.

Description of TL LUE Theory
TL LUE theory refers to stratifying a canopy into sunlit and shaded leaves groups, applying Equation 4 to calculate the GPP for each group (Chen et al., 1999).The GPP of the whole canopy is calculated as follows: where ε su and ε sh are the actual LUE of sunlit and shaded leaves, respectively; and APAR su and APAR sh are the PAR absorbed by sunlit and shaded leaves, respectively.The TL-LUE model is the first TL LUE model that was developed based on the BEPS model by separately calculating the APAR of sunlit and shaded leaves groups (Chen et al., 1999;Guan et al., 2021;He et al., 2013;Zhou et al., 2016).The RTL-LUE model was then developed based on the TL-LUE model, with the assumption that the difference in LUE between sunlit and shaded leaves is determined by the received intensity of incident radiation.In this way, the same maximum LUE is assigned for sunlit and shaded leaves, making it possible to have one maximum LUE in the same way as the BL MOD17 model (Guan et al., 2021).The GPP can be calculated as follows: where ε * max is the maximum LUE of all the leaves within the canopy; f (VPD) and g (T min ) are the same as in Equation 1; f (PPFD su ) and f (PPFD sh ) are the radiation scalars f (PPFD) for sunlit and shaded leaves, which are calculated by the reciprocal function of photosynthetic photon flux density (PPFD), as follows: where b is set as a constant with the value of 1 (mol m 2 hh 1 ) (Chen et al., 1999), and only parameter a is used to control the response of f (PPFD su ) and f (PPFD sh ) to PPFD (mol m 2 hh 1 ).The PPFD for sunlit and shaded leaves can be obtained from the APAR (PARsu or PARsh) multiplied by a constant PAR-energy ratio of 4.55 mol/ MJ (Chen et al., 1999;Running & Coughlan, 1988).The APAR sh and APAR su in Equation 5 are calculated as follows (Chen et al., 1999): where α is the canopy albedo related to the vegetation type (Table S2 in Supporting Information S1); θ is the solar zenith angle; β is the mean leaf-sun angle, which is set to 60°for a canopy with a spherical leaf angle distribution, as in previous studies (Chen et al., 1999;Guan et al., 2022;He et al., 2013); C represents the multiple scattering of direct radiation (Chen et al., 1999); PAR dif ,U is the diffuse PAR under the canopy, which is calculated following Chen et al. (1999); and PAR dif and PAR dir are the diffuse and direct components of incoming PAR, respectively, which are calculated using Equations 9 and 10: where PAR is the total incoming photosynthetically active radiation, which is equal to (S 0 × 0.5); and R is the sky clearness index, which is equal to (PAR/(0.5 × S 0 × cos θ)); S 0 is the solar constant (1367 W m 2 ).The LAI su and LAI sh are the LAI of sunlit and shaded leaves and are partitioned on the basis of the LAI in Equation 11and Equation 12.They are computed as follows (Chen et al., 1999): )) ( 11) where Ω is the clumping index, for which we use the same value as previous studies (Tang et al., 2007;Zhou et al., 2016) (Table S2 in Supporting Information S1).

Description of the CTL-MOD17 Model
The MOD17 (Equation 1) and RTL-LUE (Equation 5) models use the same environmental stress parameters f (VPD) and g (T min ), with the only difference being the vegetation canopy treatment and APAR calculation.If we divide Equation 5by Equation 1, we can obtain the following relationship: where ε max , PAR, and f PAR are calculated in the same way as in the MOD17 model.Considering the radiation difference between sunlit and shaded leaves, RPPFD su and RPPFD sh are the ratio of APAR and the radiation of sunlit and shaded leaves, respectively, which are very important for MODIS GPP correction.RPPFD su and RPPFD sh are established as follows: where APAR su , APAR sh , f (PPFD su ) , and f (PPFD sh ) are calculated in the same way as in the RTL-LUE model.We bring RPPFD su and RPPFD sh into Equation 13, and then the remaining part can be denoted as C R .C R can be approximately considered as the ratio of the LUE between the BL and TL models, as well as the impacts of f PAR and the difference of LAI in its reuse.As a result, the new correction method for the MODIS GPP product based on TL-LUE theory (CTL-MOD17) can be expressed as follows:

Model Parameterization and Evaluation
To provide the full set of correction parameters for future applications of the CTL-MOD17 model, parameters a and C R in Equations 6 and 16 were calibrated by each plant function type using GPP data from 153 FLUXNET sites (Table 2).For each site, the data from the same site-year are treated as a unit, and random 20% of them are used for validation and the remaining for parameterization.The shuffled complex evolution method of the University of Arizona was employed to implement the calibration (Duan et al., 1992;Guan et al., 2022;Zhou et al., 2016), which evaluates the model performance with the agreement index (d).This is calculated as follows: where N is the total number of measurements; P i and O i represent the corrected product (GPP_CTL) and observed EC GPP, respectively; and O and P are the mean values of the observations and predictions for all the experimental data points.
Three indicators are used in this paper to evaluate the CTL-MOD17 model performance: the determination coefficient (R 2 ), the root-mean-square error (RMSE), and the mean predictive error (bias), which are calculated as follows:

8-Day Results
The accuracy of the 8-day GPP from the original MODIS GPP product (GPP_MOD) and the product from the CTL-MOD17 model (GPP_CTL) was compared with the EC GPP for all the PFTs.The results show that a higher accuracy can be found for the CTL-MOD17 model than the MOD17 model, with a higher R 2 and lower RMSE and bias (Table 3).When analyzed for individual PFTs, the improvement of the CTL-MOD17 model was more pronounced for CSH, OSH, WSA, CRO, GRA, and EBF, with an R 2 increase of more than 0.06.For the mean of all the PFTs, the R 2 increases from 0.52 for MOD17 to 0.58 for the CTL-MOD17 model, and the corresponding RMSE and bias decrease from 17.49 to 15.80, and from 2.75 to 0.98 g C m 2 8 d 1 , representing an improvement in R 2 of +0.06 and a reduction in RMSE of 1.75 g C m 2 8 d 1 and bias of 1.76 g C m 2 8 d 1 for the CTL-MOD17 model, respectively (Table 3).The results indicate that the estimation accuracy of the CTL-MOD17 model is higher than that of the MOD17 model when considering the combined contributions of sunlit and shaded leaves to GPP.Furthermore, the corrected MODIS GPP product from the CTL-MOD17 model is more similar to the EC data, overall, with a smaller systematic error.
To further explore the overall overestimation or underestimation of MODIS GPP, we compare the scatter plots for the two models with EC GPP in 11 PFTs in Figure 1.It is apparent that the slopes of the corrected result obtained with the CTL-MOD17 model are much closer to 1 than those of the MODIS GPP product, indicating a better agreement with the EC data, except for CSH.It is apparent that the correction of the MODIS GPP product reduces overestimation in CSH, OSH, and DNF over low vegetation productivity areas and underestimation of high vegetation productivity for almost all the vegetation types.The slopes of CSH ( 0.16), OSH ( 0.08), and DNF ( 0.09) are decreased, while considerable increases in the slope can be found for CRO (+0.32),GRA (+0.31),ENF (+0.29),WSA (+0.26), and DBF (+0.22) (Figure 1).After the results for all the PFTs are averaged, the slope increases from 0.59 for the MODIS GPP to 0.82 for GPP_CTL, with a significant difference of 0.23 (Figures 1i and 1i').Overall, the CTL-MOD17 model mitigates the underestimation of the MOD17 model in the high productivity range, with a more concentrated distribution around the 1:1 line for most of the vegetation types.

Monthly and Yearly Results
To further verify the performance of the CTL-MOD17 model, the 8-day data were converted to monthly sum and yearly sum composites to reduce the temporal volatility.As shown in Figures 2a and 2b, the CTL-MOD17 model also agrees better with the EC GPP data than the MOD17 model at the monthly scale, with an increase in R 2 (+0.04) and a decrease in RMSE ( 7.38 g C m 2 month 1 ).For the yearly composite, more significant improvements can be observed in the CTL-MOD17 results.The R 2 increases from 0.71 for MOD17 to 0.78 for CTL-MOD17, corresponding to an RMSE decrease from 460.50 to 339.21 ( 121.29 g C m 2 year 1 ) and increased R 2 (+0.07) (Figures 2c and 2d).It is apparent that the GPP of the yearly composite has higher R 2 and lower RMSE than the GPP of the monthly product, because the yearly composite reduces the random errors of the daily and monthly GPP caused by the variability in environmental conditions.
For the slope, the monthly composite increases from 0.61 for MOD17 to 0.83 for CTL-MOD17, and the corresponding value of the yearly composite increases from 0.56 to 0.73, representing increases of 0.22-0.17,respectively (Figure 2).These results indicate that the CTL-MOD17 model is better at capturing seasonal changes than the MOD17 model, especially in the medium and high productivity ranges.From the scatter plots, it can be observed that the main reason for the improvements is that the CTL-MOD17 model can better simulate GPP in the high productivity range, which is often underestimated in the MOD17 results.
From the above analysis, it can be concluded that the CTL-MOD17 model can improve the accuracy of GPP simulation at 8-day, monthly, and yearly time scales.For all or individual vegetation types, the GPP corrected by the CTL-MOD17 model shows increased R 2 and reduced RMSE, relative to the MODIS GPP product, when evaluated against EC data.

Seasonal Variation of the CTL-MOD17 Model
The vegetation physiology and canopy structure change with the season can influence GPP, so it was crucial to evaluate the ability of the CTL-MOD17 model to capture the seasonal variation of GPP.The seasonal variations of GPP from the CTL-MOD17 and MOD17 models were compared with the EC GPP at  representative sites across 11 vegetation types (Figure 3).It can be found that the CTL-MOD17 model matches better with the EC GPP than the MOD17 model, with higher accuracy in capturing seasonal variations in all PFTs.All the vegetation types, except for CSH, OSH, and DNF, are underestimated in the MOD17 model, especially between DOY (day of the year) 145 and 217.In contrast, the CTL-MOD17 model captures the seasonal patterns of GPP much better in this regard.Both the original MOD17 and CTL-MOD17 results overestimate GPP for DNF, which may have been caused by the LUE saturation due to sparse leaves.
In addition, it can also be found that the simulated GPP increases significantly at the peak in the CTL-MOD17 model, compared to the MOD17 model, for CRO, DBF, EBF, ENF, GRA, and MF (Figure 3).This improvement in simulating the peak GPP with the CTL-MOD17 model can mainly be attributed to increasing the contribution of shaded leaves to GPP.

Accuracy Differences in Different Seasons
Due to the difference in the density of vegetation canopy leaves and the solar angle in different seasons, the contribution of sunlit and shaded leaves should also be different in the different seasons (Zhou et al., 2016).Therefore, we compared the differences between the "green season" and the "brown season."The "green season" represents the primary growth period of vegetation in the Northern Hemisphere (April-October), and the remaining 6 months are defined as the "brown season."The timing is the opposite in the Southern Hemisphere, compared to the Northern Hemisphere (H.Liu et al., 2013;Piao et al., 2007).In this study, since the MODIS GPP product for the DNF vegetation type lacks values for the brown season, we only compared the accuracy for the DNF green season (Table 4).For the brown season, the CTL-MOD17 model performs slightly better than the MOD17 model in most vegetation types.However, all the vegetation types show a considerable improvement with the CTL-MOD17 model during the green season, with higher R 2 and lower RMSE, especially for CSH, DBF, and OSH.After averaging all the vegetation types, the R 2 improves from 0.44 for MOD17 to 0.52 for the CTL-MOD17 model, and the corresponding RMSE and bias decrease from 19.21 to 17.28 g C m 2 8 d 1 and from 2.5 to 0.33 g C m 2 8 d 1 , respectively, representing an R 2 increase of 0.08 and an RMSE and bias decrease of 2.18-2.83g C m 2 8 d 1 , respectively, during the green season.
In the brown season, the accuracy improvement of the CTL-MOD17 model (ΔR 2 = 0.03, ΔRMSE = 0.65 g C m 2 8 d 1 , Δbias = 0.46 g C m 2 8 d 1 ) is relatively low.Furthermore, the improved accuracy of CTL-MOD17 in the green season for CRO is slightly lower than that in the brown season, which is possibly due to the human impacts (irrigation, fertilization, etc.) in the green season that are not considered in either model.
To further examine the performance of the CTL-MOD17 model relative to the MOD17 model in the green and brown seasons, we compared the total GPP estimated by the MOD17 and CTL-MOD17 models with the EC GPP in the 11 PFTs through bar graphs (Figure 4).The GPP estimated by the CTL-MOD17 model for all the vegetation types is closer to the EC GPP than the MOD17 model in the green season.The total GPP for CRO, DBF, EBF, and ENF from the CTL-MOD17 model is significantly improved.In contrast, the total GPP for CSH, DNF, and OSH from the CTL-MOD17 model is lower than for MOD17 and is closer to the EC GPP in the green season (Figure 4a).For the brown season, the performance of the CTL-MOD17 model is not as good as for the green season, and the GPP of the different vegetation types varies greatly (Figure 4b).These results show that the GPP estimation accuracy is affected by the seasonal changes and further confirm that the CTL-MOD17 model performs better than the MOD17 model, especially when the vegetation leaves are dense in the green season.

The Relationship of the Bias in GPP With LAI
Previous studies have indicated that the underestimation of the MOD17 GPP product is also related to the distribution of LAI (Zhang et al., 2012).Therefore, in order to explore whether the amount of correction  made by the CTL-MOD17 model was also correlated with LAI, the differences between the EC GPP and the original MODIS GPP product (MOD-EC) and the amount of correction by the CTL-MOD17 model (CTL-EC) were calculated and are plotted with the LAI in Figure 5.It can be observed that compared with the EC GPP, the MODIS GPP product shows obvious overestimation in the low LAI range and underestimation in the high LAI range, and the underestimation tends to increase with the increasing LAI (Figure 5a).This problem can be corrected well by the CTL-MOD17 model, as the differences between the corrected and EC GPP show an almost random distribution around the zero line.Both the R 2 and the slope between the LAI and the GPP difference are approaching zero (Figure 5b), which is a significant improvement over the MOD17 product, which has 0.1-4.45 for the R 2 and the slope.Furthermore, we also plotted the relationship between the LAI and the amount of correction, as shown in Figure 5c, where the yaxis represents the amount of correction minus the original MOD17 GPP.A more significant relationship can be found, with R 2 = 0.35 and slope = 4.96, which also indicates that the overestimation of the MODIS GPP product in low LAI regions and the underestimation in high LAI regions are both well corrected.This is mainly because the BL MOD17 model only considers the average conditions of the vegetation canopy, and thus overweights the contribution of shaded leaves in low LAI regions and underweights the contribution in high LAI regions.

Comparison of the Global Distribution Patterns
In this section, we describe how the corrected global product for 2010 was used to analyze the uncertainty distribution pattern of the CTL-MOD17 GPP by comparing it with the EC GPP and TL-LUE products and further verifying the contribution of shaded leaves to GPP estimation at a global scale.

Evaluation of the Corrected GPP Product
The performance of the CTL-MOD17 model was further verified by comparing the global distribution at different spatio-temporal scales.We generated a globally corrected GPP product at a spatial resolution of 0.05°for every 8 days for the year 2010.The accuracy of the GPP of the MOD17, TL-LUE, and CTL-MOD17 models was evaluated against the EC GPP at 104 sites (including 10 vegetation types), which have valid data in the year 2010.
The estimation results of the CTL-MOD17 model perform better than the MOD17 model in all the vegetation types, showing higher R 2 and lower RMSE and bias, especially for EBF, OSH, SAV, and CRO (Table 5).When calculating the mean value of all the vegetation types, the R 2 increases from 0.52 for the MOD17 model to 0.59 for the CTL-MOD17 model, and the corresponding RMSE and bias decrease from 20.27 to 18.09 g C m 2 8 d 1 and from 7.11 to 0.34 g C m 2 8 d 1 , with increased R 2 (+0.07) and reduced RMSE ( 2.18) and bias ( 6.77).The product of the TL-LUE model was also used for comparison, and it was found that the accuracy of the TL-LUE model is also higher than that of the MODIS GPP product, but the accuracy is slightly lower than that of the CTL- MOD17 model, especially when the bias is more obvious (Table 5).Meanwhile, the scatter plots for the MOD17 (green), TL-LUE (blue), and CTL-MOD17 (red) products at the EC sites are displayed in different colors for a direct qualitative comparison (Figure 6).The slope of the CTL-MOD17 model is closer to 1, compared to that of the TL-LUE and MOD17 models.Overall, although the TL-LUE and CTL-MOD17 models both effectively alleviate the bias of the MOD17 model at high and low LAI, and especially the underestimation problem at high LAI, the CTL-MOD17 model exhibits smaller bias ( 7.17) and RMSE ( 1.24), indicating a higher accuracy (Table 5).

Comparison of Seasonal Changes
To verify the ability of the CTL-MOD17 model to capture seasonal changes, we compared the corrected product (GPP_CTL) of the CTL-MOD17 model with the TL-LUE (GPP_TL) and MODIS GPP (GPP_MOD) products at a global 8-day scale (Figure 7) (Bi & Zhou, 2022).The results reveal that the seasonal variations of the mean GPP in GPP_CTL and GPP_TL are generally consistent, and both are higher than GPP_MOD, especially in the green season (Figure 7).Furthermore, it is also apparent that GPP_CTL provides a slightly more effective solution to the underestimation problem in MODIS GPP than GPP_TL, which may be due to the consideration of radiation effects in the CTL-MOD17 model, which has an advantage in capturing seasonal variations.These results further demonstrate the stable performance of the CTL-MOD17 model with regard to seasonal variations and its effective mitigation of the bias issue in MODIS GPP during the green season.

Comparison of the GPP Global Spatial Distribution Patterns
To determine the error distribution of the global GPP products, we also compared the global spatial patterns of the GPP products estimated from the CTL-MOD17, TL-LUE, and MOD17 models (Figure 8).It can be found that all three products have similar spatial distributions, with high GPP in low latitudes and a gradual decrease toward high latitudes.However, the GPP of the CTL-MOD17 and TL-LUE models is significantly higher than that of the MOD17 model under high LAI conditions, and these models have lower GPP in low LAI regions than the MOD17 model (Figures 8a and 8b).The CTL-MOD17 model effectively corrects the low estimation of the MOD17 model in high LAI areas, especially in the evergreen broadleaf forests near the equator, and also the high estimation in low LAI areas with shrub vegetation (Figures S1 and S3 in Supporting Information S1).Compared with the MOD17 model, the spatial distribution pattern of the CTL-MOD17 correction results is closer to that of the TL-LUE model (Figures 8c and 8d).The difference between GPP_MOD and GPP_CTL is shown in Figure 9, where significant differences can be mainly observed in the Northern Hemisphere (50°N) and equatorial regions (0°).There are three peaks and two peaks of GPP along the longitude and latitude spans, respectively, which correspond to the Amazon (60°W), Congo Basin (20°E), and Indonesia (100°E).Regardless of the latitude and longitude peaks, the GPP of the CTL- MOD17 and TL-LUE models is significantly higher than that of the MOD17 model.Furthermore, Furthermore, in order to further evaluate the spatial detail characteristics, America, China, Australia, and Brazil were selected as typical regions to compare the GPP differences after correction.The results were more obvious than at a global scale that CTL-MOD17 can effectively solve the overestimation in the region with low LAI values and the underestimation within high LAI values (Figure 10).For example, the CTL-MOD17 model effectively solves the problem of bias of MOD17 products in low LAI northern and high LAI southern regions of America (Figure 10).From the RMSE difference (ΔRMSE) between the estimated GPP of the MOD17 and CTL-MOD17 models and EC GPP, the RMSE of the CTL-MOD17 model is significantly lower than that of the MOD17 model in those regions with significant differences in GPP (Figure 9c).These results further prove that the CTL-MOD17 model integrates the contributions of sunlit and shaded leaves to GPP, and effectively solves the problem of underestimation of the MOD17 model in the areas with high LAI and overestimation in low LAI areas, thus improving the accuracy of the MODIS GPP product.

Stability of the CTL-MOD17 Model
The results shown above confirm that the CTL-MOD17 model can effectively improve the accuracy of the MODIS GPP product, but the accuracy in regions without EC sites has not been verified.Therefore, we adopted a leave-one-out cross-validation approach for each vegetation type, that is, all the data of one site were not used for the parameterization and only for the validation, and all 153 sites were looped until each site was validated.In this condition, the accuracy can be regarded as the application accuracy of the CTL-MOD17 model in an area without EC sites.As shown in Table 6, we find that the CTL-MOD17 model still shows a better spatial application accuracy than MOD17.When calculating the average value of all vegetation types, R 2 improved from 0.52 for MOD17 to 0.57 for the CTL-MOD17 model, and the corresponding bias decreased from 3.38 to 0.71 g C m 2 8 d 1 .The improvements of the CTL-MOD17 model still can be observed in all 11 vegetation types, especially in CSH, EBF, and GRA.Therefore, these results further demonstrate the spatial applicability of the CTL-MOD17 model, which is important for global applications.

Significance
Although MOD17 is the most widely used GPP product, its application is constrained by the uncertainties induced by the input data, model parameters, model structure, etc.Among these sources of uncertainties, the issues regarding input data and parameters can be hard to address and are also common for other products, but the limitations caused by the BL model structure of MOD17 can be effectively corrected.This is because many of the current TL models, which feature a superior model structure, use the same environmental stressors, and the main differences are in the recalculation of the canopy APAR.Therefore, in this paper, we have proposed the first general product-based GPP correction method (CTL-MOD17), which aims to remove the biases in the MODIS product by making use of a TL model that can simulate the variable contributions of sunlit and shaded leaves to the canopy-level GPP.The results show that the CTL-MOD17 model can not only achieve a higher GPP simulation accuracy at different time scales (Figures 1 and 2), but it can also greatly reduce the underestimation and overestimation problems of the MODIS GPP product in the regions with high LAI and low LAI (Figures 5  and 8).In terms of the global spatial pattern, a great difference in the ΔRMSE between the GPP estimated by the MOD17 and CTL-MOD17 models can be found in the area around the equator and the upper mid-latitudes with high GPP values, where the results of the CTL-LUE model are more similar to the other TL products (Figure 9).Although a plethora of BL and TL LUE models have been developed to simulate GPP, they are more or less independent of each other.This study would be the first attempt to bridge the BL and TL models, and to prove the possibility of estimating high-accuracy GPP based on an existing BL product using TL theory.Many studies have shown the superiority of TL models over BL models by comparing their individual GPP simulation results, which may include uncertainties in the respective input data processes and model parameterization.The CTL-MOD17 model proves that TL theory can consistently achieve a higher accuracy than the BL MOD17 model, without  changing any MOD17 model parameters or data inputs.This consolidates the conclusion that TL models can reduce the positive and negative biases of BL GPP in high LAI and low LAI regions.
Furthermore, this study also provides a new attempt to obtain high-accuracy GPP in an economic way based on the existing products, and without additional meteorological data inputs.The environmental stressors information has already been included in the MODIS GPP product, and even though the VPD and Tmin based scalars may be insufficient, they have been widely used in various models, especially TL models (Almeida et al., 2018;Hu et al., 2022).As a result, it is unnecessary to input environmental data into the CTL-MOD17 model because the MOD17 GPP is already input (Guan et al., 2022;Pei et al., 2022;Xie & Li, 2020).In this condition, only two parameters need to be reused, that is, the LAI and PAR, to reduce the deviation of the MODIS GPP product in the calculation of APAR.Since there are many GPP products now available, it is worth exploring whether it is possible to obtain higher-accuracy GPP estimation in a more economical way based on the existing products.Extensive efforts have already been made to improve the models used to obtain GPP, starting from the basic inputs, for example, LAI, normalized difference vegetation index, enhanced vegetation index, temperature, humidity, precipitation, CO 2 concentration, et cetera (Bi & Zhou, 2022;Zhang et al., 2017;Zheng et al., 2019).In this study, an attempt has been successfully made to develop a correction model for the MODIS GPP product by converting the BL results to be similar to those of a TL model, which extends the picture for GPP estimation in a unique way, which also has the potential to be applied to GLASS GPP and other BL-based products.In order to further prove the effectiveness of the CTL-MOD17 model, the accuracies of the RTL-LUE and CTL-MOD17 models were also compared, based on the EC GPP at 153 sites.As shown in Figure S4 and Table S3 in Supporting Information S1, it can be observed that, over all the vegetation types, CTL-MOD17 exhibits a slightly higher accuracy than the RTL-LUE model, overall, with lower bias ( 0.53) and RMSE ( 0.5), and also across the different vegetation types.As a result, the CTL-MOD17 model does not need to use all the parameters required by a TL model, but it can achieve better accuracy than the RTL-LUE model, which is a finding that will further advance global terrestrial carbon cycle research.

Uncertainties
Although the correction results of the CTL-MOD17 obtained in this study appear satisfactory, only the bias in the structure of the MOD17 model is considered, and there are also many other sources of uncertainties in the MOD17 GPP product that need to be further addressed in the future (Guan et al., 2022;Pei et al., 2022;Zhang et al., 2012).These uncertainties can come from the unreliable input data, the scale differences of the multi-source data, the complications of background reflectance, the limitations of the environmental stressors, and the vegetation typespecific parameters of the used model (Baldocchi & Penuelas, 2019;Pei et al., 2022;Serbin et al., 2013;Wagle et al., 2016).With regard to the input data, many studies have already indicated that the flaws in the MODIS LAI and MERRA-2 PAR data are a key source of error, and especially the errors in the LAI data would carry over into the MOD17 product (Fang et al., 2012;Liu et al., 2018;Serbin et al., 2013).In this study, in order to avoid the impacts of the input data, exactly the same LAI and PAR inputs were used as for the MOD17 product, and thus all the improvements in GPP accuracy can be attributed to the change in model structure.However, it will also be important to explore whether another LAI product can achieve better results, and therefore we also used the GLASS LAI instead of the MODIS LAI to drive the CTL-MOD17 model.The accuracy of the correction results can be found in Table S4 in Supporting Information S1, which indicates similar overall results to those obtained based on the MODIS LAI, but with differences in the different vegetation types.This demonstrates that there are impacts caused by the input LAI data, but these are difficult to correct, except by recalculating the GPP obtained by the TL LUE models using other LAI inputs, which was not the purpose of this research.The condition would be the same for the impact of the PAR data, and thus the same LAI and PAR input data as MOD17 GPP were used in this study.Furthermore, the evaluation of the model and GPP products can also be impacted by the heterogeneity of the ecosystems and environmental conditions, which can induce huge scale differences between the used meteorological and remote sensing data, and it will be important to analyze the impacts of land cover and environmental heterogeneity on the evaluation results (Chen et al., 2013;Xie et al., 2021;Zhang et al., 2012).The spatial resolution of the PAR, MODIS LAI, and vegetation type data used in this study was not only mismatched with the footprints of the EC flux towers, but was also not entirely consistent with the surrounding vegetation types, potentially also impacting the accuracy of the GPP estimation (Pei et al., 2022;Xie et al., 2021;Zhou et al., 2016).
Furthermore, the vegetation type-specific parameters and insufficient environmental stressors in the MOD17 model may also bring uncertainties to the MOD17 GPP product, which will also need further study.The calibration of C R can be regarded as solving some of the uncertainties induced by the used lookup table of maximum LUE in the MOD17 model, because C R contains the difference between the BL and TL maximum LUE.However, the parameters in the environmental scalars (i.e., VPDmax, VPDmin, T max , T min ) were not considered in this study, which are the default values inherited from the BIOME-BGC model (Running et al., 2000), and were not parameterized using field measurements.Research has suggested that these default parameters introduce a certain level of uncertainty, and optimizing these parameters could further improve the precision of the GPP simulation (Huang et al., 2021).The VPD and Tmin parameters in the environmental scalars in the MOD17 model have also been shown to be insufficient to describe the environmental stressors.Some studies have also indicated that it could be possible to further improve the accuracy of the GPP simulation by fully considering the environmental factors, such as soil nutrients, stand age, and CO 2 (Lacroix et al., 2022;Xie et al., 2020;Zheng et al., 2018), which could also be a direction for future studies.Using the continuous CI data rather than the vegetation type-specific constants would also be a direction to further improve the correction accuracy, although it would need additional data inputs, because our results indicate that the effective LAI (LAIe = LAI*CI) can show a higher correlation with the corrected error of MOD17 GPP (Figure S3 in Supporting Information S1).This study has proved the feasibility of obtaining high-accuracy GPP by correcting the BL-based MOD17 product using TL theory.There are also many other BL-based GPP products, such as GLASS GPP, EC GPP, and VPM GPP (Zhang et al., 2017), and it will be interesting to explore the performance when applying the CTL-MOD17 model to these other GPP products.Due to these GPP products being driven by different LAI and PAR data, and the fact that the environmental scalars of the applied models are also different, comparing the performance of the correction results for these different GPP products could also help to answer the above question, to a certain extent, regarding the impacts of the input data and environmental scalars in the model.Finally, this study was aimed at improving the traditional LUE models.In recent years, SIF and ML models have emerged as new technologies capable of more accurately estimating GPP and accounting for environmental uncertainties (Bai et al., 2022;H. Li et al., 2022).Therefore, coupling LUE models, SIF, and ML methods could well be an important direction for accurately simulating GPP in the future.

Conclusion
In this study, we developed a simple model named the CTL-MOD17 to improve the accuracy of the widely used MODIS GPP product by bridging BL and TL models.This new model solves the issue of MODIS GPP underestimation for dense vegetation by introducing the capability of a TL model that can properly consider the contribution of shaded leaves to the total canopy GPP.The performance of the new model was assessed based on 153 global EC sites.The main conclusions can be drawn as follows: 1.The CTL-MOD17 model can significantly improve the accuracy of the original MODIS GPP product.The CTL-MOD17 model (R 2 = 0.78, RMSE = 339.21g C m 2 year 1 ) obtained significantly higher R 2 and lower RMSE than the original MOD17 product (R 2 = 0.71, RMSE = 460.50g C m 2 year 1 ), when compared with EC GPP at 8-day, monthly, and yearly time scales, for all 11 vegetation types.Furthermore, the CTL-MOD17 model (ΔR 2 = 0.05, ΔRMSE = 0.38 g C m 2 8 d 1 ) was found to have good stability and robustness when validated against the same EC data set using a leave-one-out cross-validation method.2. A greater improvement of the CTL-MOD17 model over the MODIS GPP product was found in the green season than in the brown season.The results of the CTL-MOD17 model can also better capture the seasonal variation in GPP and greatly alleviate the underestimation issue for the MODIS GPP product in the green season, thus achieving a much higher accuracy in the green season.This effectively minimizes the errors caused by seasonal variations and enhances the reliability of the MODIS GPP product.3. The MODIS GPP product has a serious underestimation problem in high LAI areas and a serious overestimation problem in low LAI areas.Through the ability of the CTL-MOD17 model to consider the contribution of shaded leaves to the canopy-level GPP, the issues of underestimation and overestimation of the MODIS GPP product in high and low LAI areas can be well addressed, resulting in values that are much closer to the EC GPP and other TL model products.These results demonstrate the importance of separate simulation of sunlit and shaded leaf GPP to obtain canopy-level GPP.

Figure 1 .
Figure 1.The accuracy of the 8-day GPP estimated by the MOD17 model (GPP_MOD) and CTL-MOD17 model (GPP _CTL), compared with the corresponding tower measurements (EC GPP) for all the vegetation types in the validation set.

Figure 2 .
Figure 2. The relationship between eddy covariance Gross primary productivity (GPP) and the corrected GPP results from the monthly sum (a) of the Moderate Resolution Imaging Spectroradiometer GPP product and (b) the CTL-MOD17 model, and for the yearly sum (c) and (d) composites.

Figure 3 .
Figure 3. Seasonal variations in the Gross primary productivity (GPP) simulated by the MOD17 and CTL-MOD17 models for different vegetation types, compared with the eddy covariance GPP at representative sites.The title of each subplot is named as "vegetation types/site year/site name."DOY means the day of the year.

Figure 4 .
Figure 4.The mean Gross primary productivity (GPP) between the eddy covariance GPP and the simulation results from the CTL-MOD17 and MOD17 models in each plant function type.(a) Green season.(b) Brown season.

Figure 5 .
Figure 5.The relationship between the Gross primary productivity (GPP) difference (ΔGPP) and leaf area index.(a) MOD-EC is the difference between the eddy covariance (EC) GPP and the Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product, (b) CTL-EC is the difference between the EC GPP and the CTL-LUE results, and (c) CTL-MOD is the CTL-MOD17 correction results minus the corresponding MODIS GPP product.The units of the RMSE and bias are g C m 2 8d 1 .

Figure 8 .
Figure 8. Maximum annual composite of the leaf area index (LAI) of vegetation and a comparison of the spatial distribution characteristics of the 8-day product composited into a yearly global-scale Gross primary productivity (GPP) product for 2010.(a) LAI.(b) Moderate Resolution Imaging Spectroradiometer GPP product.(c) GPP product from the CTL-MOD17 model.(d) TL-LUE product.Figures (b) and (c) are the same as the legend of (d).

Figure 9 .
Figure 9. Spatial distribution of the differences in the annual summed Gross primary productivity (GPP) correction for 2010.(a) The graphs show the differences in longitude from the MOD17, CTL-MOD17, and TL-LUE models.(b) Graph for latitude.(c) ΔGPP means GPP of the CTL-MOD17 model minus the corresponding value of the MOD17 model; ΔRMSE refers to the RMSE of the estimated results of the CTL-MOD17 model and the eddy covariance GPP of the site minus the corresponding value of MOD17, so the negative points with larger values mean better correction results.

Figure 10 .
Figure 10.The distribution characteristics of the 2010 annual summed Gross primary productivity (GPP) of the CTL-MOD17 model (CTL) and MOD17 model (MOD) products for America, China, Australia, and Brazil, and the relationship between the GPP difference of the two models (ΔGPP) and maximum leaf area index.The unit is g C m 2 year 1 .

Table 1
Number of EC Sites and Site-Years for the Different Vegetation Types in the IGBP Classification System

Table 2
The C R , a, and Agreement Index for the Different Vegetation Types a The abbreviations for the different vegetation types are the same as those in Table1.

Table 3
Statistics of the Gross Primary Productivity Simulation Accuracy for the Different Vegetation TypesThe abbreviations for the different vegetation types are the same as those in Table1.bTheunits of the RMSE and bias are g C m 2 8d 1 . a

Table 4
GPP_MOD Accuracy Statistics for the Different Vegetation Types in the Green Season and Brown Season a The abbreviations for the different vegetation types are the same as those in Table1.bTheunits of the RMSE and bias are g C m 2 8d 1 .

Table 5
Comparison of the Accuracy of the Products From the MOD17, TL-LUE, and CTL-MOD17 Models With 104 EC SitesThe abbreviations for the different vegetation types are the same as those in Table1.bTheunits of the RMSE and bias are g C m 2 8d 1 . a

Table 6
Statistics for the Stability of the CTL-MOD17 ModelThe abbreviations for the different vegetation types are the same as those in Table1.bTheunits of the RMSE and bias are g C m 2 8d 1 . a