Integration and Comparison of Multiple Two-Leaf Light Use Efficiency Models Across Global Flux Sites

Accurate estimation of gross primary productivity (GPP) from the regional to global scale is essential in modeling carbon cycle processes. The recently-developed two-leaf light use efficiency (TL-LUE) model and its revised versions based on different concepts have significantly improved the underlying mechanisms between model assumptions and photosynthetic processing. Yet few studies have compared the advantages of the various two-leaf LUE models for their practical applications. Here, an integrated model referred to as a three-parameter radiation-constrained mountain TL-LUE (RMTL3-LUE) is proposed by combining the radiation scalar of the [radiation-constrained TL-LUE model] and the topographic parameters of the [mountainous TL-LUE model]. In this way, the importance of light intensity and topography on vegetation photosynthesis is integrated. Our calibration and validation of RMTL3-LUE were carried out for 11 ecosystems with in situ eddy covariance measurements around the globe. This indicates that the model can effectively improve the GPP estimates compared with its predecessors. At the landscape scale, RMTL3-LUE can also realistically quantify topographic effects on photosynthesis, with topographic sensitivities of decreasing (increasing) with the slope on the unshaded (shaded) terrain. Furthermore, RMTL3-LUE displays an asymmetric sensitivity to PAR variability, with a low sensitivity to PAR compared with other models under high PAR conditions and a similar sensitivity to PAR in low PARs. Altogether, it is clear that the integration of the merits of multiple TL-LUE models can further improve the photosynthetic processes for various conditions amid more challenges in constructing more complex models.


I. INTRODUCTION
G ROSS primary productivity (GPP) of terrestrial ecosystems, defined as the amount of carbon fixed by plants through photosynthesis for a given unit of time and space, is the main process of the carbon cycle between land and atmosphere systems. Therefore, accurate modeling of GPP at different spatiotemporal scales is a prerequisite to understanding the terrestrial carbon cycle [1]. It is an active focal point of global change studies, especially in the decades since eddy covariance (EC) flux towers [2] and satellite remotely sensed data [3] became available.
The EC technique can provide a near real-time measure of the net exchange of CO 2 (i.e., net ecosystem exchange, NEE) between the land surface and the atmosphere with high accuracy, through which GPP can be directly estimated by differentiating ecosystem respiration and NEE [4]. Although its footprint only covers an area of a few hundred meters to kilometers, a network of global EC towers has provided key ancillary datasets for the calibration and validation of GPP models in different ecosystems [5], [7], [8]. In recent decades, many models have been developed to simulate GPP at scales from local to global, with substantial improvements in model structure and quality of forcing inputs. These models are generally classified as process models [9], [10], [11] and light use efficiency (LUE) models [12], [13], [14], [15], [16], [17]. Because process models have a complex structure and require many parameters for vegetation, soil, and climate, LUE models have been preferred for some applications. This is primarily due to the successful use of the remotely sensed vegetation index to characterize absorbed photosynthetically active radiation (APAR) [18], [19], [20], [21], [22], [23]. Earlier remote sensing LUE models simply assumed the vegetation canopy to be a large extended leaf (big-leaf concept) based on the leaf-level Monteith function (i.e., big-leaf LUE model) [24]. They would indirectly estimate GPP by multiplying APAR, maximum LUE, and environmental stress scalars at the ecosystem level. Unfortunately, these models (e.g., CASA, CFlux, VPM, MOD17) ignore the differences in APAR and maximum LUE between sunlit and shaded leaves of the canopy [6], [25], [26]. For any ecosystem, sunlit leaves receive This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ both direct and diffuse photosynthetically active radiation (PAR) and have low LUE due to light saturation. In contrast, shaded leaves absorb only diffuse PAR and have high LUE (two-leaf concept) [20], [23], [27]. Treating sunlit and shaded leaves the same in big-leaf LUE models thus led to substantial biases in GPP simulations, with underestimates on cloudy days and overestimates on sunny days [20], [28], [29]. In contrast, the two-leaf LUE model (e.g., TL-LUE) estimates GPP by APAR and LUE independently for the sunlit and shaded leaves. It has been proven to be more accurate than big-leaf LUE models (e.g., MOD17 model) in six plant functional types (PFTs) in China [20] and nine PFTs across the globe [27]. As a key parameter of TL-LUE, the LUE for sunlit and shaded leaves is calculated with an independent maximum LUE value for the sunlit (ε max-sun ) and shaded leaves (ε max-shd ), and then downscaling the ε max-sun and ε max-shd by identical environmental scalars of temperature and atmospheric water vapor deficit (VPD) (see methods). This means that specification of ε max-sun and ε max-shd is crucial for the LUE estimates of sunlit and shaded leaves.
A recent study by Guan et al. [30] found that the divergence in LUE between sunlit and shaded leaves is closely related to their radiation interception abilities rather than the differences in maximum LUE. This is because the physiological traits of a leaf are almost the same irrespectively of "sunlit" or "shaded," i.e., ε max-sun equals ε max-shd . Based on this assumption, Guan et al. [30] proposed a radiation-constrained TL-LUE (i.e., RTL-LUE) model by introducing a radiation scalar for the differences in LUE between sunlit and shaded leaves, which resulted in more accurate estimates of GPP across 169 globally distributed EC sites [30]. Unfortunately, both TL-LUE and RTL-LUE neglect the topographic effects on vegetation photosynthesis. However, the topographic attributes (elevation, slope, aspect, etc.) are known to change vegetation interceptive capability of radiation by altering the geometric relationship between the canopy and sunlight [31]. For example, Chen et al. [10] showed that the net primary productivity was overestimated by 5% in boreal mountain forests when topographic effects were neglected. Xie and Li [32] proposed a two-leaf LUE model for mountainous areas (MTL-LUE) to address this shortcoming by including topographic effects on direct and diffuse radiation and sunlit canopy area. This model not only improves the GPP estimates for mountainous forests but also presents the spatial distribution of the effects of rugged terrain on photosynthesis [32], [33]. The MTL-LUE, however, neglects radiation-related differences in LUE between sunlit and shaded leaves addressed by RTL-LUE, likely leading to an overestimation of actual LUE and GPP. Additionally, the MTL-LUE model has been validated at the site-or watershed-scale in only a few areas (e.g., forests in the mountains of central Italy and southwestern China) and has not been extensively assessed at EC sites globally.
Previous studies have provided important information on methodological advances in the two-leaf LUE models in global EC sites and terrestrial ecosystems. However, to our knowledge, several methodological challenges remain when integrating three two-leaf LUE models (i.e., TL-LUE, RTL-LUE, and MTL-LUE). First, both TL-and RTL-LUE models neglect the topographic effects on APAR in GPP estimates. However, both emphasize the differences in ε max (TL-LUE) and radiation intensity (RTL-LUE) between sunlit and shaded leaves. To address this, we first propose a radiation-constrained TL-LUE model for mountain areas [hereafter RMTL2-LUE] by incorporating APAR estimates of MTL-LUE with consideration of topographic effects into RTL-LUE (see methods). The inspiration for this treatment is that the RTL-LUE simulates GPP more accurately than TL-LUE by quantifying different radiation scalars for sunlit and shaded leaves. Second, we propose a radiation-constrained TL-LUE model for mountain areas that considers a different ε max for sunlit and shaded leaves (hereafter RMTL3-LUE). For example, the RMTL3-LUE incorporates APAR estimates in MTL-LUE and differing radiation constraints for the LUE of sunlit/shaded leaves of RTL-LUE into the TL-LUE (see methods). This integrated model involves the physiologically based process of photosynthesis (i.e., different ε max values for sunlit and shaded leaves) and radiation constraints on photosynthesis (i.e., different radiation intensities for sunlit and shaded leaves) under complex terrain conditions (i.e., different estimates of APAR for flat and mountain areas). We expect that RMTL3-LUE could more precisely capture the photosynthesis processes of the canopy. However, comprehensive comparisons among these models and validation against ground measurements of GPP (e.g., EC-based observations) are lacking, particularly across various PFTs and terrain conditions. Finally, the sensitivity of these models to crucial forcing inputs (PAR, LAI, and clumping index) and terrain variability is needed to examine which forcing input is most associated with the outputs of the models under diverse topography conditions. Sensitivity analysis can also identify the possible effects of biases in the gridded products of PAR, LAI, and clumping index data on spatial distributions of GPP in mountainous areas. Here we leverage 658 site-year high-quality EC data (143 sites, 11 PFTs globally) to calibrate the above two integrated models at the PFT scale and compare them with their predecessors (TL-, RTL-, and MTL-LUE). We focus on which of these models is the most effective, especially in different terrain conditions. Our premise is that new integrated models could effectively explore the merits of earlier two-leaf LUE models and not only improve the GPP estimates at EC sites but also reveal the topographic effects on spatial distributions of GPP, particularly for the mountain ecosystems at large spatial scales.

II. DATA AND METHODS
A. Data and Sources 1) Flux Data: Flux data recorded at 143 EC sites of the global FLUXNET 2015 database were used to calibrate and validate five types of models, including our two integrated models and their three predecessors (TL-, RTL-, and MTL-LUE).
The FLUXNET 2015 provides data on GPP, air temperature, precipitation, PAR, etc., at multiple time scales (30 min to annual) [34]. The high-quality flux and micrometeorological data were first selected based on the criteria, with (1) record of daily GPP, downward shortwave radiation, VPD, and air temperature, (2) daily data quality flag >0.8, (3) proportion of high-quality data for the whole year >80%, and (4) availability of remotely sensed leaf area index (LAI) and digital elevation model (DEM) data in a rectangular area of 1000 m × 1000 m centered on the EC site (hereafter carbon footprint area). The 8-day flux data for each site was then estimated by averaging the daily data with valid values of ≥5 days during consecutive 8-day records. This was to match the 8-day LAI data. As a result, the 143 EC sites (658 site-years) were selected from all 212 sites to cover the 11 PFTs globally: deciduous broadleaf forest (DBF, 19 sites), evergreen broadleaf forest (EBF, 14 sites), evergreen needleleaf forest (ENF, 28 sites), mixed forest (MF, nine sites), grass (GRA, 27 sites), crop (CRO, 15 sites), closed shrub (CSH, two sites), open shrub (OSH, seven sites), wetlands (WET, ten sites), savannas (SAV, six sites), and woody savannas (WSA, six sites) ( Fig. 1 and Table S3). Finally, the 658 site-years of flux data were randomly divided into 497 and 161 site-years for calibration and validation of all models (75.5% and 24.5% of all site-years), respectively.
2) Leaf Area Index: The Terra MODIS LAI product MOD15A2H, with a spatial resolution of 500 m and 8-day interval, was used to drive the LUE models in the study. Spatially averaged LAI time series data were extracted over each EC site's carbon footprint area after removing cloud-contaminated pixels using quality control flags. Then, a locally adjusted cubic spline capping filter was used to smooth the LAI data further to reduce the influences of nonvegetation effects such as clouds, snow, and soil [35].
3) Digital Elevation Model: The NASADEM, a newly released product generated from the original Shuttle Radar Topography Mission data with 30 m spatial resolution, provides key topographic parameters for driving the topographic-related LUE models (MTL-LUE and our integrated models). The topographic variables needed in the models include slope, aspect, and sky view factors to provide the proportion of sunlit and shaded leaves in the canopy and the distribution of solar radiation (i.e., direct and diffuse radiation).

B. Model Improvement Protocols 1) Model Description and Integration:
The model improvement protocols in this study were inspired by the methodological advances from the TL-LUE to RTL-LUE and MTL-LUE. Specifically, the TL-LUE estimates GPP by setting a different maximum LUE and APAR for sunlit and shaded canopies (1) shown at the bottom of the next page, and assuming different physiological traits of sunlit and shaded canopies. In contrast, RTL-LUE emphasizes radiation constraints rather than physiological mechanisms on LUE for sunlit and shaded canopies and employs different radiation scalars to downscale the maximum LUE of two canopies (2) shown at the bottom of the next page. However, these two types of models are generally parameterized for flat areas without consideration of topographic effects on solar radiation received by the canopies. The MTL-LUE model estimates GPP using a terrain-adjusted APAR for sunlit and shaded canopies in mountain areas, with a similar parameterization of maximum LUE in TL-LUE [20] (3) shown at the bottom of the next page. Here, we propose two integrated models for EC sites and landscape within the carbon footprint area [31] by taking the merits of the previous three models. In the first protocol, we propose a model (RMTL2-LUE) that integrates the terrain-corrected APAR of MTL-LUE into the RTL-LUE (4) shown at the bottom of this page. As its name indicates, this model estimates different APARs for sunlit and shaded canopies in mountainous vegetation and has two PFT-specific parameters to be optimized, including the light intensity coefficient [a] and ε max . In the second protocol, we propose a model (RMTL3-LUE) by integrating the terrain-corrected APAR of MTL-LUE and the radiation scalar of RTL-LUE into the TL-LUE (5) shown at the bottom of this page. Three PFT-specific parameters are to be optimized, including the light intensity coefficient [a], and the maximum LUE for sunlit and shaded canopies [ε msu and ε msh ], where the light intensity coefficients were inherited from optimized values of RMTL2-LUE. where ε max-sun , ε max-shd , and ε max are the maximum LUE of sunlit, shaded, and the entire canopy, respectively (unit: g C MJ −1 ). APAR sun-ft (APAR sun-mt ) and APAR shd-ft (APAR shd-mt ) are the absorbed PAR by sunlit and shaded leaves in flat (rugged) terrain (unit: MJ·m −2 ·d −1 ). The specific information and associated parameters optimization are provided in Supplementary Methods S1 and S2. According to Chen et al. [9] and Guan et al. [30], the radiation scalars [f(PPFD sun ) and f(PPFD shd )] were derived from the algorithms of stomatal conductance in response to light intensity in the BEPS and BIOME-BGC models. These are the reciprocal equations of the photosynthetic photon flux density (PPFD) of the sunlit and shaded leaves.
where a and b are the coefficients that determine the relationship between light intensity and LUE; b is set as a constant (1 mol m −2 hh −1 ), and a is a parameter that controls the response of LUE to PPFD. The PPFD of sunlit and shaded leaves (PPFD sun and PPFD shd , mol m −2 hh −1 ) was derived from the product of the PAR absorbed by the two sets of leaves with a constant PAR-energy ratio (4.55 mol/MJ), respectively. f(VPD) and g(T a ) represent water and temperature stress scalars that are used to downscale maximum LUE (8) and (9) where VPD and T a are day-time averaged vapor pressure deficit and daily minimum temperature in 8-day steps, respectively. The subscripts max and min are the corresponding maximum and minimum values of VPD and T a and are determined by the PFTs following Guan et al. [30].

2) Model Parameterization and Validation:
The maximum LUEs [ε max-sun , ε max-shd, and ε max ] and the radiation scalars [f(PPFD sun ) and f(PPFD shd )] are the key parameters in the five LUE models to be optimized using the EC GPP for each PFT.
In this study, the shuffled complex evolution-University of Arizona method (SCE-UA) [36] was used to optimize biomespecific parameters for each model, including the maximum LUE [ε max-sun , ε max-shd, and ε max ] and a (for f(PPFD) estimate) in (1)-(5), respectively. Specifically, SCE-UA optimization was carried out at 497 site-years and evaluated with agreement index (d), (10), which ranges between 0 (no agreement) and 1 (complete agreement). Site-year data with d <0.5 was excluded from the modeling calibration to eliminate the effects caused by abnormal or incomplete seasonal changes in EC GPP.
where n is the total number of observations; P k and O k represent the kth model predicted, and EC site observed GPP values, respectively, andŌ represents the mean value of all observations.
The statistics of optimized parameters of the five LUE models for 11 PFTs are shown in Table S4. Three statistical indexes, including coefficient of determination (R 2 ), root-mean-square error (RMSE), and mean predictive error (Bias), were applied to assess the accuracy of the 8-day GPP simulated by the five models [30] using validation data at 161 site-years.

C. Simulations Across Heterogeneous Spaces
To evaluate model performance in mountainous areas, we selected six mountainous sites (also representing six PFTs, Fig. 1) with an average slope >10°within the carbon footprint areas to simulate 8-day GPP (Figs. 5 and 6). The simulations were made at a 30-m spatial resolution to match the spatial resolution of the topography data. The forcing inputs for each model over the footprint area [T a and VPD in 8-day steps] were obtained from the Mountain Microclimate Simulation Model 4.3 (MT-CLIM 4.3) by inputting topography and flux tower micrometeorological data [37], [38], [39]. Direct and diffuse PARs above the canopies within the tower footprint were calculated using the method described in Supplementary Method S1. We uniformly used mean LAI data within the footprint areas to model landscape scale GPP for each pixel in the footprint areas. This is because the 30 m resolution LAI data for each pixel over the footprint area were unavailable, and carbon footprint areas were found to be relatively homogeneous by our visual inspection of Google Earth for each mountain site.

III. RESULTS
A. Multiscale Comparison Among the Models 1) EC Site Scale: After calibration of the models using 497 site-years data, the estimated EC GPP data at 161 site-years were used to validate each two-leaf LUE model for the 11 PFTs at an 8-day interval (Fig. 2). Overall, all five LUE models reasonably captured the variations of GPP at EC sites for 11 PFTs, with R 2 of >0.50 for all PFTs except CSH (0.213-0.315). Compared with the TL-LUE, both RTL-LUE and MTL-LUE showed better performance for most PFTs. For example, the smallest RMSE values are observed for RTL-LUE, with ranges from 5.868 g C m −2 8d −1 for OSH to 21.379 g C m −2 8d −1 for CRO (Fig. 2). The phenomena of over-or under-estimates for low and high GPP values in MTL-LUE were significantly reduced, with the fitted line more closely located at the 1:1 line than TL-LUE for most PFTs (Fig. 2). In terms of our integrated models, RMTL3-LUE outperformed RMTL2-LUE for most PFTs, with an R 2 of 0.289-0.853 and 0.315-0.872. An exception occurred in croplands where R 2 is 0.642 and 0.652 for RMTL3-LUE and RMTL2-LUE, respectively. Because RMTL3-LUE effectively integrated the merits of its predecessors, RMTL3-LUE produced a more accurate GPP prediction than TL-LUE and RTL-LUE. However, it was slightly inferior to MTL-LUE in solving the phenomena of over-and under-estimation of GPP in the low and high GPP ranges. It is worth noting that RMTL3-LUE produced higher R 2 values than predecessors for DBF, EBF, ENF, CSH, OSH, and WET, but similar for other PFTs except cropland. When all PFTs were considered together (Fig. 2), RMTL3-LUE showed the best performance, with the highest R 2 being 0.761. Additionally, RMTL3-LUE simulated GPP and flux tower GPP exhibited high levels of consistency in both the high-and low-value ranges. This was similar to MTL-LUE, superior to RTL-LUE and RMTL2-LUE, but inferior to TL-LUE.
Statistics at independent EC sites showed that more than 60% of the sites produced an R 2 of >0.7 and an RMSE of <12 g C m −2 8d −1 , and nearly 80% of the sites exhibited |Bias| <10 g C m −2 8d −1 for all models (Table I). This suggested that all models were universally applicable in various regions of the world (Fig. 3). In particular, RMTL3-LUE produced the highest R 2 values and the lowest RMSE and |Bias| values at 47.6%, 46.3%, and 36.6% of the validation sites, which is considerably better than the other four models (Table II).  To evaluate the performances of the five models in rugged terrain, we divided EC sites into two groups by the mean slope in the footprint area. The first group represents the flat terrain where the mean slope is 0°-5°(64 sites, 112 site-years), while the second group represents rugged terrain where the mean slope is >5°(18 sites, 39 site-years). Interestingly, we found that the five models produced similar results for the first group, with R 2 , RMSE, and Bias ranging from 0.765-0.770, 17.742-15.507 g C m −2 8d −1 and −2.808-−1.753 g C m −2 8d −1 . In addition, the regression lines between the observed GPP and modeled GPP by TL-LUE, MTL-LUE, and RMTL3-LUE were very close to the 1:1 line with a k of >0.8 and a b of <7.6 g C m −2 8d −1 (Fig. 4). In rugged terrain, however, the five models' modeling accuracy was reduced to varying degrees, with a noticeable reduction of the R 2 and k values for TL-LUE (from 0.765 to 0.718 and from 0.836 to 0.758) and the k value for RMTL2-LUE (from 0.783 to 0.749). The RTL-LUE, MTL-LUE, and RMTL3-LUE showed more consistency with EC GPP in rugged terrains, most probably due to the successful integration of radiation/topography mechanisms in these models. In summary, RMTL3-LUE showed similar or better performance than other models regardless of PFTs, sites, and terrain features. However, RMTL2-LUE did not although it also considers topography effects on GPP simulations. Thus, we excluded RMTL2-LUE in subsequent analysis.
2) Landscape Scale: Clear spatial differences were observed in total GPP (i.e., the sum of sunlit and shaded GPP) between modeled by [TL-LUE, RTL-LUE] and [MTL-LUE, RMTL3-LUE] (Fig. 5, S2, and S4). Specifically, the GPP modeled by TL-LUE and RTL-LUE were higher in low elevations, while GPP simulated by MTL-LUE and RMTL3-LUE were higher in unshaded terrain (Figs. 5 and 6). For example, the regional averaged total GPP modeled by RMTL3-LUE at the CN-Din site was 51.544 g C m −2 8d −1 in unshaded terrains, which is greater than that in shaded terrains (49.019 g C m −2 8d −1 ) (Fig. 6). As    6. Spatial distribution of the annual mean GPP simulated by four LUE models in varying aspects (denoted by radius) and slopes (denoted by angle) at six EC sites (i.e., six PFTs). The black dashed line indicates the annual mean solar azimuth (As); the red/blue number represents the spatial average for annual mean GPP in unshaded terrain [As-90°<aspect< = As+90°]/shaded terrain [As+90°<aspect< = 360°or 0°<aspect< = As-90°]. expected, both MTL-LUE and RMTL3-LUE showed greater spatial variabilities than TL-LUE and RTL-LUE, with the standard deviations of modeled total GPP ranging between 0.132 and 2.202 g C m −2 8d −1 for MTL-LUE and 0.101 and 1.941 g C m −2 8d −1 for RMTL3-LUE (Figs. 5 and 6). This suggests that topographic corrected models can reveal the spatial variations of GPP. In terms of the spatial differences in GPP between sunlit and shaded canopies (GPP sun and GPP shd ), the nontopographic-corrected models (TL-LUE and RTL-LUE) showed similar spatial patterns of GPP sun and GPP shd with total GPP for all selected landscapes, with a tendency of increasing GPP sun and GPP shd with altitudes (Figs. S2-S5, Figs. 5 and 6). However, both MTL-LUE and RMTL3-LUE exhibited different spatial GPP sun and GPP shd , with higher GPP sun in unshaded terrain than in shaded terrain, which is consistent with the spatial pattern of total GPP and higher GPP shd in shaded terrain than in unshaded terrain (Figs. S2-S5, Figs. 5 and 6).

1) Sensitivity of GPP Simulated by Models to Key Input
Variables: We optimized maximum LUE and light intensity coefficients using EC data globally and emphasized the advantages of integrating the merits of previously developed two-leaf models. However, the variations of modeled GPP could be associated with three key inputting variables, including LAI, PAR, and clumping index. Consequently, we further analyzed the sensitivity of the model output (GPP) to the above three variables for each PFT by changing these variables relative to the original simulation data from −30% to 30% in steps of 5%. The sensitivity of all models to PAR tended to decrease sharply with increasing PAR (from 0% to 30%) and to increase slightly with decreasing PAR (from 0% to −30%) for most PFTs, showing an asymmetric distribution of the logarithmic function (Fig. 7). For example, a 30% increase in PAR resulted in a 12.1% overestimation of GPP, and a 30% decrease in PAR resulted in a 19.8% underestimation of GPP in MTL-LUE for EBF. It is noteworthy that the sensitivity of RMTL3-LUE to PAR is the lowest for all PFTs, with GPP variations ranging from −22.0% to +16.2%, which is lower than that of other models (−30.9% to +42.0%; Fig. 7). Interestingly, there are almost no differences in the sensitivity of all models to LAI for individual PFT. However, the LAI sensitivity varied greatly by PFTs, with high sensitivity occurring in CRO, OSH, SAV, and WSA (−27.4% to +25.4%) and low sensitivity in the other PFTs (−20% to +20%) (Fig. 8).
The sensitivity of the four models to the clumping index showed similar patterns to the LAI (i.e., individual PFT). There was a slight difference in sensitivity between models, with differences of GPP variability of <11.7%, while for the independent model, there was a larger difference between PFTs, with GPP variability exceeding 20% in OSH, SAV, and WSA, but −16-11% for other PFTs (Fig. 9).

2) Sensitivity of the GPP Simulated by Models to Topographic Variations:
In addition to three inputting variables, the topography is another factor causing variations of GPP modeled by the topography-corrected models (MTL-LUE and RMTL3-LUE). By artificially changing the slope from 0 to 30°in 1°s teps for unshaded and shaded terrain separately, we showed that the GPP simulated by both RMTL3-LUE and MTL-LUE exhibited higher sensitivity to slopes on shaded terrain than unshaded terrain and that sensitivity increased with an increasing slope on shaded terrain but decreased with the slope on unshaded terrain (Fig. 10). Compared with MTL-LUE, the GPP simulated by RMTL3-LUE showed a lower sensitivity to the slope on unshaded terrain and a higher sensitivity on shaded terrain. Considering that the input variables in different terrain may play a decisive role in the final GPP variations by the topography-corrected LUE model, we further found that the sensitivity of RMTL3-LUE and MTL-LUE to changes in key input variables with varying terrain. We found that the sensitivity of the two models to PAR and LAI increased with slope increase on unshaded terrain and decreased with the slope on shaded terrain (Figs. S6 and S7). In contrast, the sensitivity of the two models to the clumping index showed opposite changes with the slope on the unshaded or shaded terrain (Figs. S8). For the same topographic conditions (specific slope and aspect), RMTL3-LUE presented lower sensitivity to PAR, higher sensitivity to LAI, and similar sensitivity to the clumping index compared with MTL-LUE (Figs. S6-S8).

A. Model Improvement
Using EC and micrometeorological data in 11 PFTs from 143 flux sites globally, and the modeling concept of two-leaf LUE models, we integrated the merits of three previously proposed two-leaf LUE models [TL-LUE, RTL-LUE, and MTL-LUE] into two individual models. In the first integration, we optimized the PFT-dependent maximum LUE and light intensity coefficients in the entire canopy without consideration of differences in ε max between sunlit and shaded canopies. In the second integration, we optimized the different ε max for sunlit and shaded canopies, except for the light intensity coefficient optimized by our first integration. In addition, all two integrations also estimated the APAR for sunlit and shaded canopies separately by considering the effect of rugged terrain on APAR. Comparing the model estimates using validation data in 86 EC sites showed that RMTL3-LUE (second integration) provides a noticeably better model performance than RMTL2-LUE (first integration) for all PFTs. This suggests that the independent setting of ε max for sunlit and shaded canopies is necessary for improving the prediction of GPP. Notably, the RMTL3-LUE produced similar or more accurate estimates of GPP than three original two-leaf models, especially in the forest (DBF, EBF, and ENF) and shrubland (CSH and OSH) ecosystems (Fig. 2), indicating broader applicability of the newly integrated model. This improved estimate of GPP in RMTL3-LUE is most probably associated with the full consideration of differences between both ε max and light intensity for sunlit and shaded canopies, and the effects of topographic dependency on canopy solar radiation reception. For example, He et al. [20] suggested that a separate specification of ε max for sunlit and shaded canopy can increase the accuracy of GPP estimation compared with the big-leaf model MOD17 that use only one ε max for the whole canopy. Guan et al. [30] also increased the accuracy of GPP by separately parameterizing the radiation constraint (light intensity) on ε max for the sunlit and shaded leaves. More importantly, relevant field experiments and mechanistic studies further support the above model mechanisms. Boulard et al. [40], Yoshimoto et al. [41], Miller et al. [42], and Chen et al. [43] observed that different micrometeorological conditions at the sunlit and shaded canopy interact with canopy structure to influence canopy physiology (enzyme activity, diffusion rate, etc.). This physiological difference between sunlit and shaded canopies may lead to the necessity of setting ε max for these two groups of leaves separately. Hirose and Werger [44] also proposed that the plant canopy adapted to the radiation environment by redistributing nitrogen to maximize photosynthesis throughout the canopy. They also proposed that photosynthesis is stronger in the upper layers with higher nitrogen content than in the lower layers with lower nitrogen content, supporting the rationality of different radiation constraints to constrain the ε max in sunlit and shaded canopies in our model. The improvements of RMTL3-LUE can also be observed at mountain sites compared with the MTL-LUE model proposed by Xie and Li [32] (Fig. 4).
Despite using the same topographic parameters, RMTL3-LUE captured the GPP of mountain vegetation more accurately than MTL-LUE, which may be attributed to the integration of merit of RTL-LUE (radiative constraints) to the RMTL3-LUE. From sensitivity analysis, it can also be observed that RMTL3-LUE is less sensitive to PAR (Fig. 7) and significantly less variable with increasing slope than MTL-LUE (Fig. S6), indicating the attenuation of uncertainty in mountain PAR simulation to GPP simulation in RMTL3-LUE model. The gridded PAR products were reported to overestimate the actual PAR [46]. In this respect, the low sensitivity of RMTL3-LUE to PAR, especially for high PAR values, may improve the accuracy of regional-scale GPP simulations. At the landscape scale, RMTL3-LUE produced more reasonable spatial distributions of GPP, GPP sun, and GPP shd than TL-LUE and RTL-LUE, consistently to MTL-LUE (Figs. 5 and 6, Figs. S2-S5) [32], [33]. This divergent spatial pattern of GPP modeled by the RMTL3and MTL-LUE model is attributed to the spatial simulation of the model-driven variables of LAI and PAR considering rugged terrains. By comparing the spatial patterns of LAI and PAR from RMTL3-LUE with results from Xie and Li [32], Wang et al. [45], and Yu et al. [46] at six mountain sites, we found that the LAI and PAR estimated for mountainous areas by RMTL3-LUE in conditions of terrain-tilting and adjacent terrain shading were generally consistent with the results of previous studies. Specifically, sunlit LAI for RMTL3-LUE increases with an increasing slope on unshaded terrains and decreases with an increasing slope on shaded terrains, while the opposite occurs for shaded LAI (Figs. S9 and S10). It is observed that the direct PAR for RMTL3-LUE is higher on unshaded terrains than on shaded terrains, where direct PAR increases with an increasing slope on unshaded terrains and decreases with an increasing slope on the shaded terrains. Conversely, diffuse PAR decreases with increasing slope regardless of the aspect variations (Figs. S11 and S12). However, the TL-LUE and RTL-LUE models, without consideration of topographic parameters, failed to model the spatial variations in LAI and PAR (Figs. S9-S12) and the final GPP simulations for different terrains (Figs. 5 and 6). The minor spatial variation of GPP with elevation presented by them was induced by the MT-CLIM 4.3 algorithm's elevation correction factor for temperature, independent of the internal parameters of the TL-LUE and RTL-LUE models. At the site scale, RMTL3-LUE effectively captured GPP at nearly 70% of globally distributed EC sites (R 2 > 0.7, RMSE < 12 g C m −2 8d −1 , and |Bias| < 10 g C m −2 8d −1 at 68.3%, 69.5%, and 80.5% of validated sites in this study, Table I). RMTL3-LUE outperformed other two-leaf models at nearly 50% of EC sites (RMTL3-LUE had the highest R 2 , lowest RMSE, and |Bias| at 47.6%, 46.3%, and 36.6% of validated sites, Table II), demonstrating the robust performance of RMTL3-LUE in capturing terrestrial GPP in different geographical regions of the globe. Although improvements in GPP simulation at both site and landscape scales were observed in the RMTL3-LUE model, several key concerns remain regarding the application of the RMTL3-LUE model in the future. The first concern is the increased complexity of the model structure and, consequently, the requirements of larger input data and computation cost for large regional scale applications than the predecessors (TL-, RTL-, and MTL-LUE). Due to the integration of three previously proposed models, three PFT-specific parameters, including ε max-sun , ε max-shd , and a, needed to be optimized. Second, we limited our study to the flux sites and less than 20 mountain sites due to the difficulties in parameterizing models in areas without carbon flux and micrometeorological observations. An establishment of global scale ε max-sun , ε max-shd, and a will become necessary to improve or apply the model per pixel of remote sensing data already available globally. Therefore, the trade-off in application between RMTL3-LUE and its predecessors depends mainly on the requirements of study using these models. Suppose the regional or local terrain heterogeneity has a decisive effect on GPP variation, and there is a need for precise GPP simulation. In that case, we strongly encourage employing the RMTL3-LUE model to explore the effects of micrometeorology on carbon fluxes. If the studies are on the continental or global scales that may allow neglecting some small-scale heterogeneity, the selection of the models highly depend on the availability of flux data, topography data, and the study requirement.

B. Model Parameterization
The light intensity coefficient [a] and the maximum LUE [ε max , ε max-sun , and ε max-shd ] for sunlit and shaded canopies are pivotal parameters optimized independently before the application of the two-leaf LUE models. Overall, the a, ε max , ε max-sun , and ε max-shd values optimized here varied with model structures, site-year, and PFT (Table S4). In agreement with the studies by He et al. [20], Zhou et al. [27], and Guan et al. [30], the mean maximum LUE decreased with an order of ε max-shd > ε max > ε max-sun for each model and PFT, with a mean value of 2.44±0.36, 1.74±0.90 and 1.12±0.50 g C MJ −1 , respectively. Similarly, the coefficient of light intensity [a] in RTL-LUE and our integrated models [RMTL2-LUE, RMTL3-LUE] are also consistent for PFTs, with values ranging from 0.18 to 1.76 and 0.22 to 1.96, respectively. However, there are some variations in these optimized parameters for an individual model in different site-years, suggesting the importance of ecosystem representation and high-quality site-years EC data for flux towers. For instance, the standard deviations of a, ε max-sun and ε max-shd for RMTL3-LUE were 0.20-1.07, 0.34-1.88 g C MJ −1 and 0.48-1.68 g C MJ −1 among the PFTs, and the coefficients of variation were 20.6%-115.4%, 29.7%-69.6% and 13.6%-51.6% (Table  S4). This suggests that the high-quality site-years flux data are essential for model parametrization, and setting fixed values for specific PFT may increase the uncertainty of the ultimate GPP simulation. Therefore, the accuracy of GPP prediction could be improved enormously if these parameters were reoptimized well using EC observations, especially in precise small-scale studies. If EC observations are unavailable, a relatively robust GPP simulation (Fig. 2) can be obtained by directly using the PFT-specific parameters optimized in this study (Table S4).
Accurate simulation of PAR in mountainous areas is key to achieving accurate estimates of mountain GPP at regional scales. In nontopographic corrected models [TL-and RTL-LUE], the PAR was divided into direct and diffuse components (PAR dir and PAR dif ) based on the empirical relationships between PAR and the sky clearness index (SI). In contrast, the topographic corrected models [MTL-and RMTL3-LUE] simulated PAR dir and PAR dif by SI-based transmittance of direct and diffuse radiation (K b and K d ) and topographic auxiliary inputs (Eqs. S7 and S8) [32]. We recalibrated K b and K d in this study at 42 EC sites globally (Table S1). We showed that both PAR dir and PAR dif simulated by topographic corrected models using K b and K d are highly in agreement with those observed in validation sites, with [k = 0.997, R 2 = 0.896] for PAR dir and [k = 1.048, R 2 = 0.791] for PAR dif , respectively. In comparison, the nontopographic corrected models underestimate PAR dir and overestimate PAR dif , with [k = 0.622, R 2 = 0.769] for PAR dir and [k = 1.696, R 2 = 0.795] for PAR dif (Fig. S1). In summary, establishing the relationship between corrected SI and K b (K d ) for the global multiecosystems is imperative to achieve accurate GPP estimates for RMTL3-LUE.

V. CONCLUSION
This study developed an integrated two-leaf LUE model (RMTL3-LUE) by integrating the radiation scalars of RTL-LUE and the topographic parameters of MTL-LUE into TL-LUE and considering the differences in their effects on sunlit and shaded canopies. Comparisons of RMTL3-LUE with its predecessors showed that the newly proposed model increased GPP predictions for most plant functional types and EC sites. This highlighted that integration of the merits of different models could be important, especially for areas that need precise estimates of GPP for exploring micrometeorological effects on carbon cycles. However, the model structure may become complex in the new model. For complex terrain conditions, RMTL3-LUE is more effective in resisting the attenuation effect of tilted and obscured terrain on model performance than other two-leaf models. It realistically characterizes the spatial patterns of GPP in mountainous landscapes. Sensitivity analysis of the model inputs shows that RMTL3-LUE is significantly less sensitive to PAR than other models, which could serve to resist the unfavorable effects of the current overestimated grid radiation data inputs on GPP simulations. We do not observe significant intermodel divergence regarding sensitivity to LAI or clumping index in this study. Detailed studies using more precise LAI and clumping index data with a finer spatial resolution are expected in the future. Considering the full optimization of three key PFT-specific parameters, including a, ε max-sun , and ε max-shd , we speculate that the new integrated model involves more biophysical mechanisms of plant photosynthesis. For example, canopy physiological attributes and radiation constraints on photosynthesis for sunlit and shaded canopies, and differences in APAR for unshaded and shaded terrain. Our model optimization work for RMTL3-LUE provides key equations related to mountain PAR allocation and PFT-specific parameters over 11 ecosystems globally. It is expected that RMTL3-LUE is capable of more correctly quantifying terrestrial productivity from landscape to watershed level in different geographical areas globally.