Distinguishing the main climatic drivers to the variability of gross primary productivity at global FLUXNET sites

Climate exerts both short-term and long-term impacts on the ecosystem carbon assimilation. However, the main climatic drivers for the variability of gross primary productivity (GPP) remain unclear across various timescales and vegetation types. Here, we combine the state-of-the-art machine learning algorithms with a well-established explanatory method to explore the impacts of climatic factors on long-term GPP variability at global FLUXNET sites across four timescales and six plant functional types. Results show that diffuse shortwave radiation (SWdif) dominates GPP variability at the sub-daily (half-hourly to three hourly) timescales especially for the tree species, and acts as the secondary contributor after air temperature at the daily or longer timescales. Attribution analyses further showed that the main effects of SWdif are much higher than their interactive effects with other climatic factors in regulating the GPP variability. By identifying the main climatic drivers, this study improves the understanding of the climate-driven GPP variability and provides important implications for the future projection of ecosystem carbon assimilation under global climate change.


Introduction
Terrestrial ecosystems absorbed over 30% of anthropogenic carbon dioxide (CO 2 ) emissions in the past several decades, thus mitigating global climate warming (Friedlingstein et al 2022).The strength of such carbon sink is initially associated with the gross primary productivity (GPP) that determines the total amount of carbon assimilation by plant photosynthesis (Fernández-Martínez et al 2017, Ma et al 2017).Observational and modeling studies have revealed large spatiotemporal variabilities in GPP (Yue et al 2015, Zhang and Ye 2021), which are caused by both biotic and abiotic factors including climatic conditions (Sun et al 2018), nutrient supply (von Sperber et al 2017), plant species (Moore et al 2017), and physiological functions (Chapin et al 2002).Among those factors, climatic conditions act complex roles through both short-term meteorological variations and long-term climate change (Piao et al 2013, Niu et al 2017), leading to large uncertainties in the estimation of present-day sinks (Sitch et al 2015) and the projection of future changes in ecosystem carbon fluxes (Tian et al 2021).
Previous estimates of the climatic impacts on GPP variability showed varied conclusions for different timescales (van Dijk et al 2005, Yuan et al 2007).At the sub-daily to daily timescales, climatic factors with strong diurnal variations may drive GPP variability (Han et al 2020).At the longer timescales, the impact of meteorology can be masked by environmental factors with slow but accumulated effects on the long-term period (Wu et al 2017).For example, Ishtiaq and Abdul-Aziz (2015) found that shortwave radiation (SW) had 5-8 times stronger impacts than temperature on the half-hourly carbon fluxes at eight deciduous forest sites in the U.S. In contrast, Cao et al (2017) showed that temperature explained 83%-93% of monthly GPP variability at an alpine wetland site in China.Even for the same site, different meteorological factors may exert dominant impacts at multiple timescales.For example, Mäkelä et al (2006) found that SW explained up to 95% of the sub-daily GPP variability at a boreal forest site, but temperature instead of radiation dominated GPP variability at the longer timescales.Meanwhile, water availability could also module plant photosynthesis through root moisture absorption (Green et al 2019) or stomatal open and closure (Yuan et al 2019).At the short-term timescales, atmospheric vapor pressure deficit (VPD) usually controls GPP variability as an important factor via instantaneous modulations of plant stoma and photosynthesis (Fu et al 2022).For example, Wang et al (2022b) found that VPD contributes the highest importance to GPP variability at the growing seasons using datadriven models.At the longer timescales, soil moisture instead of VPD dominates plant carbon assimilations by altering total water availability in plant roots (Green et al 2019).Furthermore, changes in phenology, length of growing season, and managements may also influence GPP variability across various timescales.
Apart from timescales, the impacts of climatic factors on GPP vary across various plant functional types (PFTs).For instance, Zhou et al (2011) found that SW is the most critical factor affecting halfhourly carbon fluxes at two poplar forest sites in China.For the same half-hourly timescale, however, Lee et al (2015) found that temperature instead of radiation controlled plant photosynthesis at two wetland sites in China.For radiation, canopy structure is usually the key factor affecting light availability and the consequent leaf photosynthesis (Chen et al 2022).For temperature, the heat tolerance and leaf age are instead the more important factors determining the GPP variability.Such differences explain why the GPP variability of tropical forests, which have strong heat tolerance, are more sensitive to radiation and water constraints than temperature (Madani et al 2020).Furthermore, leaf nitrogen content, specific leaf weight, and chemical properties of the vegetation are also key factors influencing relationships between climatic changes and terrestrial carbon uptake (Valentini et al 1999, Reichstein et al 2014, Migliavacca et al 2021).Considering these findings are derived based on limited observations using varied methods, it remains unclear how the climatic factors influence global GPP variability (main and interactive effects) across different timescales and PFTs.
In addition, diffuse radiation also plays an important role to terrestrial carbon uptake due to radiative disturbances of aerosols and cloud (Hutchison and Hicks 1985, Jarvis et al 1985, Hollinger et al 1994, Freedman et al 2001, Gu et al 2002, Alton et al 2007).Observations have shown that diffuse radiation can help increase the light use efficiency of canopy and therefore promote ecosystem GPP (Farquhar and Roderick 2003, Gu et al 2003, Cirino et al 2014, Yang et al 2019, Zhou et al 2022).Nevertheless, the impact of diffuse radiation on GPP variability is not fully explored due to the lack of diffuse radiation datasets with wide spatial coverage and high temporal frequency (Zhou et al 2021c).Among the limited explorations, Han et al (2020) analyzed flux data at a forest site using the path analysis method and found that diffuse radiation modulates daily variability of GPP to a larger extent than other climatic factors such as direct radiation and temperature.Similarly, Gui et al (2021) explored the diffuse fertilization effect on GPP at nine sites in China and showed the dominant role of diffuse radiation in driving the GPP variability at half-hourly timescales.These studies suggest that diffuse radiation may act as the lead driver to GPP under certain circumstances.
In this study, we explore the linkages between climatic factors and the observed long-term GPP variability at the global FLUXNET sites by combining the state-of-the-art machine learning algorithms and a newly-developed explanatory method.We aim to (1) examine how climatic factors regulate GPP variability across various timescales and PFTs, and (2) identify the main and interactive effects of climatic factors on modulating GPP variability.We specifically focus on exploring the role of diffuse radiation neglected in previous studies, taking advantage of a recently developed dataset of hourly diffuse radiation at global flux tower sites (Zhou et al 2021c).We also compare the contributions of the main and interactive effects so as to understand how the associations between climate factors may jointly affect the GPP variability.

Observations
FLUXNET (http://fluxnet.fluxdata.org/) is a worldwide network measuring ecosystem fluxes based on eddy covariance methods.It provides long-term observations of carbon and water fluxes, and simultaneous meteorological measurements at half-hourly timescales (Pastorello et al 2020).In this study, we select 52 sites with at least ten year observations of carbon fluxes from FLUXNET2015 product, including six main PFTs (table S1 and figure S1).We use the carbon fluxes of GPP (µmol CO 2 m −2 s −1 ), concentrations of carbon dioxide (CO 2 , ppm), and meteorological variables of air temperature (TA, • C), SW (W m −2 ), soil water content (SWC, %), VPD (hPa), surface pressure (P, hPa), wind speed (WS, m s −1 ), and precipitation (PRE, mm hr −1 ).To reduce uncertainties of net ecosystem exchange (NEE)-based methods, GPP data was used as the mean of estimations between the daytime (GPP_DT_VUT_MEAN) and nighttime (GPP_NT_VUT_MEAN) partitioning methods.Since the deep SWC was observed at very few sites (Zhou et al 2019), we use the derived SWC data at top soil layer in the FLUXENT2015 product integrated by the marginal distribution sampling method.Continuous observations of diffuse fraction are not available at most FLUXNET sites.Here, we use the site-level diffuse fraction derived using an artificial neural network (ANN) approach combined with observations and reanalyses datasets (Zhou et al 2021c), and calculate diffuse shortwave radiation (SWdif, W m −2 ) and direct shortwave radiation (SWdir, W m −2 ) using derived diffuse fraction and observed SW.Validations show that the ANN approach can reproduce observed temporal variations of diffuse fraction with higher correlations and lower biases than 11 other empirical methods.The earlier applications also showed its good capability in exploring the impacts of diffuse radiation on ecosystem carbon (Zhou et al 2021b) and water fluxes (Wang et al 2022a).The diurnal and seasonal variations of climatic factors averaged across all the global FLUXNET sites are shown in figures 2 and S3.We analyze the simultaneous carbon and meteorological variables only if over 80% of daytime (SW >0) records are available.In this study, we define peak growing seasons as June-July-August for sites in northern hemisphere and December-January-February for sites in southern hemisphere.

Extreme gradient boosting (XGB) model
Gradient boosting machine (GBM) is a representative tree-based model that reduce biases between training and observed datasets using gradient descent algorithms step by step (Newton and Wernisch 2017).However, GBM often causes overfitting due to the effect of hyperparameters and quality of training datasets (Agapitos et al 2017).As an updated algorithm, XGB overcomes the overfitting by using parallel processing and regularized boosting (Chen and Guestrin 2016).Different from traditional linear regression methods, XGB algorithm uses decision trees to construct predictive models based on inputoutput relationships.Decision trees are often by nature immune to co-linearity between these input variables, due to selection of the optimal feature.For example, once there are two features with 99% of correlation, decision tree will choose only one of them as statistic feature within the processes of deciding upon a split.Recently, the XGB method was used to address problems in atmospheric environment (Park et al 2021) and terrestrial ecosystems (Wang et al 2022b) and showed better potentials than traditional methods.Here, we apply the XGB model to reproduce site-level GPP based on eight climatic drivers including TA, diffuse radiation (SWdif), direct radiation (SWdir), VPD, CO 2 , WS, air pressure (P) and PRE.For each site, four GPP-XGB models are trained with climatic factors averaged at half-hourly, three-hourly, daily and monthly timescales, respectively.Evaluations show that these GPP-XGB models reproduce observed GPP with correlation coefficients higher than 0.85 across various timescales for 90% FLUXNET sites (table S2).

Shapley additive explanation (SHAP) method
SHAP method is an explainability technique to quantify the contributions of each input variable to the predicted values based on the theory of cooperative games (Lundberg and Lee 2017).This method directly regulates the effects of each feature on model losses or gains using the differences between predicted and expected values (Fryer et al 2021).Recently, the XGB model was applied as a foundation for the SHAP-based interpretations (Lundberg et al 2020), and was used to address problems in public health (Nohara et al 2022), socioeconomic sciences (Mao et al 2021), and ecosystems (Wang et al 2022b).In this study, we combine the SHAP method with the XGB model to quantify the contributions of individual climatic factors and their interactive effects to GPP variability across different timescales and PFTs at global FLUXNET sites.
For each sample x with N climatic features, the SHAP method is used to calculate SHAP values of each input variable as follows (Aas et al 2021): where f and g are predicted and explanation models, respectively.Simplified input x ′ is used to replace x to match explanation models with the original predicted model.φ 0 is the mean value of predicted GPP, and φ i x ′ i is the SHAP values for climatic feature i in sample x.The sum of SHAP values for all features and the mean predicted GPP is equal to the predictions by the XGB model.Negative/positive SHAP values indicate inhibiting/promoting effects on GPP caused by climatic features in sample x.Different from the importance index used in other methods, such SHAP values are calculated for each feature in each sample, providing the opportunity to explore the contributions of climatic factors to diurnal or seasonal variabilities of GPP.
For each feature, SHAP values can be divided into main and interactive effects based on shapley interaction index (Grabisch andRoubens 1999, Fujimoto et al 2006) as follows: where

Impacts of climatic factors on GPP variability for various timescales
The SHAP values of each factor in the prediction of GPP are quantified at various timescales (figure 1 S5).Not all the SHAP values of climatic factors are positive (figure 1(b)).High SWdif, TA, and SWdir usually induce positive SHAP values.SHAP values for SWC first increase along with low SWC, and then stabilize to zeros with higher SWC (figure S6), indicating non-linear controlling effects of soil water availability.Meanwhile, slight increases in VPD cause positive SHAP values, which may be related to interactions between plant transpiration and photosynthesis.However, the higher VPD often causes the larger GPP losses (negative SHAP values) owing to stomatal closures (figure S6).Such GPP responses to climatic factors are consistent with previous empirical responses derived using regressions (Brandao et al 2013) and simulations (Yuan et al 2019).CO 2 and other factors have no significant impacts on GPP variability as shown by the narrow range of SHAP values.For different timescales, SWdif and TA always act as the top two factors regulating GPP variability (figure 1(c)).The importance of SWdif decreases from half-hourly to monthly timescales, while the contributions by TA increase gradually and become the dominant factor at the daily to monthly timescales.Meanwhile, SWC accounts for stronger impacts on GPP variability than that by SWdir at the daily and monthly timescales, indicating higher importance of soil water availability for the long-term variability of plant photosynthesis.Furthermore, contributions of other climatic factors remain stable from half-hourly to daily timescales, except for PRE which shows larger contributions at longer time scales.For the monthly timescale, our results show effects of climatic factors are largely dampened especially for SWdif, SWdir, and VPD, indicating that climate-induced GPP variability is smoothed during the time accumulation.

Impacts of climatic factors on GPP variability for various PFTs
We explore the SHAP values of climatic factors in the prediction of GPP for various PFTs (figure 2), and our results show controlling effects of SWdif on tree species are higher than that of non-tree species.At the half-hourly timescale, SWdif dominates GPP variability for all PFTs except deciduous broadleaf forest (DBF) (figure 2(a)).For tree species, SWdif contributes 25.6%-39.2% to the GPP variability, higher than the level of 22.9%-28.2%for non-tree species.As the secondary driver, TA makes contributions of 5.7%-29.7% to GPP variability across six PFTs with the minimum for evergreen broadleaf forest (EBF) and the maximum for DBF.For the rest factors, SWdir, SWC and VPD make stable contributions of 12.7%-29.8%,6.3%-20.8%,and 5%-10.7%,respectively, and other factors show limited impacts on GPP with low SHAP values.At the daily timescale, TA on average contributes 33.5% to the mean absolute SHAP values across all PFTs (figure 2(b)), followed by the contributions of 21.5% by SWdif, 12% by SWC and 9% by SWdir.In general, TA dominates the GPP variability of most PFTs except for EBF, Crop and Shrub; for two former PFTs, the contributions of 8.9% and 15.2% by TA is lower than that of 26% and 24.3% by SWdif.For Shrub, however, TA contribute mean absolute SHAP values by 17.9%, lower than 24.3% caused by SWC.From half-hourly to daily timescales, our results show the contributions by SWdif and SWdir reduce by 8.1% and 6.5% while those by TA increase by 11.2%, with the largest shifts during the growing season for most PFTs (figure S7).As a comparison, other factors show no significant changes among various PFTs between the two timescales.With the longer timescales, contributions of SWdif and SWdir averagely decrease up to 15.9% and 5.9% especially for tree species (EBF, evergreen needleleaf forest (ENF) and DBF) at the monthly scales (figure S8).Instead, TA increased by 44.7% across all PFTs, contributed by 59.6% for ENF and 64.4% for DBF.Other climatic factors show limited changes among most PFTs from half-hourly to monthly scales.S2.

The main and interactive effects of climatic factors
We further separate the contributions of climatic factors into main and interactive effects (figure 3), and our results show main effects of climatic drivers especially for SWdif and TA are stronger than interactive effects of other factors among various timescales.For the main effects, SWdif shows the highest SHAP contribution of 21.4% at the half-hourly timescale, followed by the contributions of 16.9% by TA, 10.1% by SWdir, 6.5% by SWC, 6.1% by VPD, and 4.8% by CO 2 (figure 3(a)).These main effects of individual factors are much larger than their interactive effects, which show the highest value of 3.3% for SWdif-TA associations and the secondary value of 2.1% for TA-VPD associations.Similarly, the main effects instead of the interactive effects of climatic factors contribute to most of SHAP values at the daily timescale (figure 3(b)), but with the largest SHAP contribution of 26.1% by TA and the secondary contribution of 16% by SWdif.With the shift of timescales from half-hourly to daily, the SHAP contributions to GPP variability by climatic factors show reductions of 5.3% for SWdif, 5% for SWdir, and 1% for VPD, but increase 9.2% for TA and 0.95% for SWC (figure 3(c)).At the monthly timescale, the contributions by the main effects of TA increase up to 34.2% while SWdif and SWdir contribute only 13% and 4% (figure S9).For interactive effects, our results show SWC-and PRE-related associations gradually enhance from the half-hourly to longer timescales (figure 3), indicating accumulated effects of soil water in the processes of long-term carbon assimilation.However, the GPP contributions by VPDrelated associations decreased continually from halfhourly to monthly timescales, indicating that atmospheric water conditions caused stronger effects on GPP variability at the short-term instead of long-term timescales.
We examine the main and interactive effects of climatic factors on GPP variability for various PFTs (figure 4).At the half-hourly timescale, the main effect of SWdif results in the highest mean absolute SHAP values of 16.5%-30.1%for six PFTs, followed by the main effects of TA, SWdir, SWC and VPD (figure 4(a)).The main effect of SWdif plays the dominant role in regulating GPP variability for EBF, ENF, Shrub and Crop sites, but shows comparable contributions to the main effect of TA for DBF and Grass sites.At the daily timescale, the main effect of TA dominates GPP variability with SHAP values of 11.1%-33.6%for all PFTs but EBF (figure 4(b)).
For EBF, although the SHAP value of SWdif is largely dampened by −11.6% at the daily timescale, its magnitude is 3.5 times that of TA at the same timescale.For all PFTs, the shift of GPP timescale from half-hourly to daily results in increases of 0.9%-12.5% in the SHAP values for TA but reductions of 2.5%-11.6%for SWdif and 3.2%-11.3%for SWdir (figure 4(c)).Meanwhile, the main effects of SWC increased by 2.9%-4.6%from half-hourly to daily timescales at EBF, DBF and Shrub sites, but showed limited contributions for other PFTs.With the longer timescales, the main effects by TA increase continually up to 34.6%-52.9%for ENF, DBF and Grass sites at the monthly timescales (figure S10).For EBF, SWdir (18.1%) and SWC (12%) play more important roles than SWdif (12%) and TA (4.4%) in driving GPP variability.For shrubland, SWdif and SWC are still the dominant GPP drivers at the monthly timescales.Compared to main effects, the interactions among climatic factors show limited impacts on GPP variability for various PFTs from half-hourly to monthly timescales.

Comparisons with the literature
We collected a total of 30 researches through the literature review (table S3) and summarized the dominant climatic drivers of carbon fluxes from these studies (figure 5).It was found that solar radiation usually dominated the GPP variability at hourly to daily timescales, while temperature dominates GPP variability at daily to longer timescales.This finding indicates ecosystem GPP are more sensitive to light availability for the daily or shorter timescales.For plant, solar radiation directly determined the availability of light energy used for instantaneous photosynthesis (Medlyn et al 2005), thus the strong diurnal cycle of solar radiation (from ∼0 in the morning to >500 W m −2 at the noontime) results in the variation of photosynthesis from 0 to the maximum within the day.This range of photosynthetic rate is much larger than that induced by temperature, which usually varies within 10 • C in the daytime and causes a variation of photosynthesis less than 20% (Farquhar et al 1980).For the longer timescales, however, impacts of solar radiation on ecosystem GPP decreases gradually owing to the smoothed diurnal variability in the radiative fluxes.Instead, TA or SWC, PRE became dominant GPP drivers with increased seasonal variations and the accumulated warming/cooling or wetting/drying effects on plant phenology (Fernandez-Martinez et al 2014, Ma et al 2019).
Meanwhile, the impacts of climatic factors on GPP also vary across various PFTs (figure 5), indicating importance of canopy structures to plant carbon uptake.We found that SW always dominates plant photosynthesis for tree species with complex canopy structures, while temperature becomes more important for GPP variability of non-tree species.However, the net impact of climatic factors on GPP variability sometimes depends on the competition between the effects of timescale and PFT.For example, Zhou et al (2009) estimated GPP drivers at the half-hourly timescale for a wetland site and found that temperature rather than radiation controls the GPP variability, indicating that the effects of PFT are more important than those of timescale in northeastern China.On the contrary, Han et al (2020) found that SW instead  S3.
of temperature dominated daily GPP variability at a subtropical coniferous forest site.
Here, we present comprehensive assessments of climatic effects on GPP variability by combining long-term observed carbon fluxes of global FLUXNET sites with the well-established XGB-SHAP framework.Our results show that shortwave radiation especially its diffuse portion (SWdif) dominates GPP variability over 30% at half-hourly/hourly timescales (figure 1) with more significant impacts for tree species (figure 2).Meanwhile, our derived responses showed TA instead of radiation as the main GPP drivers at the longer timescales (figure 1), and these findings are consistent with previous estimates (figure 5).Moreover, our results further revealed that main effects of climatic factors determine GPP variability at the half-hourly and daily timescales (figures 3 and 4), much higher than interactive effects of climatic factors such as SWdif-TA, SWdir-TA and TA-VPD.Similarly, Gui et al (2021) also found that SWdif controls GPP variability by direct instead of indirect effects over a few sites based on the path analysis method at the half-hourly timescale.
Compared to the literature, our study provides the comprehensive assessment of GPP drivers on the global scale with following advances: (1) previous studies examined GPP drivers at the single or a few sites with limited time spans (usually 1-5 years), while this study explored half-hourly data for at least ten years over 52 FLUXNET sites and revealed more universal relationships between GPP and climatic factors.(2) The XGB-SHAP framework successfully reproduced the observed GPP variability with correlation coefficients higher than 0.85 for most of FLUXNET sites, and as a result provided a more robust estimate of the contributions of individual climatic factors to GPP variability.(3) With the data-driven method, the linkage between GPP and climatic factors was explored at different temporal scales (hourly to monthly) and cycles (diurnal and seasonal) across varied vegetation types (six PFTs) for both the main and interactive effects, which were neglected in previous studies.

Importance of diffuse radiation on GPP variability
Our study emphasized the large contributions of SWdif to GPP variability across various timescales and PFTs.The SWdir promotes the photosynthesis of sunlit leaves to a saturation state that the further increase of direct light no longer enhances carbon assimilation.As a comparison, SWdif regulates the photosynthesis of shading leaves whose photosynthesis is usually in the light-limited regime (Mercado et al 2009, Wang et al 2018).As a result, the ecosystem GPP may be more sensitive to the changes of SWdif within the diurnal cycle due to the stronger photosynthetic responses (Zhou et al 2021c).Within the diurnal cycles, SWdif has higher impacts (>30%) on GPP in the sunrise, noontime, and sunset periods (figure S3).At the sunrise and sunset, high diffuse fraction drives GPP variability because SWdif is usually larger than SWdir (Zhou et al 2021a).At the noontime, SWdif with higher light use efficiency dominates GPP variability through the regulation of shading leaves, while the higher SWdir shows limited impacts due to the light-saturation state of sunlit leaves (Alton et al 2007).For the daily timescales, SWdif on average contributes 20% of GPP variability with relatively higher importance in the summer (high diffuse light use efficiency) and winter (high diffuse radiative fraction) (figures S2 and S4).
The difference of PFTs is also an important factor affecting GPP drivers especially for radiation.Our results showed that the importance of SWdif was higher for tree species than for non-tree species (figure 2(a)).Similarly, Zhou et al (2021c) and Chen et al (2022) found that SWdif caused stronger GPP responses for tree species than non-tree species.Such difference in PFT responses may be related to the complex canopy structure with more light-limited leaves for tree species (Majasalmi and Rautiainen 2020).As a result, trees have the larger fraction of shaded leaves relying on the SWdif for carbon assimilation than the non-tree species.We also found that TA caused secondary contributions to GPP variability for all PFTs except EBF.The exceptional response to temperature for EBF is likely due to the inhibition effect of high ambient temperature in the tropics on GPP variability (Tian et al 2021).Our analyses reveal that solar radiation, SWC and VPD make higher contributions to GPP variability of EBF than TA (figure 2(b)), consistent with the findings by Aguilos et al (2018) and suggesting that the carbon uptake of tropical forest is more limited by sunlight and water than heat.Furthermore, our study emphasizes importance of canopy structure as a key trait parameter distinguishing the SWdif impacts on GPP between forest and non-forest species.Other trait parameters (e.g.leaf age, specific leaf weight) may also influence the climatic effects on GPP variability across various timescales and deserve further investigations with more data availability in the trait features.
By isolating main and interactive effects, our results show SWdif present positive responses of main SHAP values (figure S11), indicating the dominant role of main effects relative to interactive effects across various timescales and PFTs.For interactions, SWdif show significant responses of interactive SHAP values to TA and VPD, but no significant interactions with SWC (figure S11).Among hourly timescales, the interaction effects between SWdif and TA, VPD display diverse patterns.SWdif-TA interactions always present positive SHAP values with lower (higher) SWdif and TA, and vice versa.Meanwhile, SWdif-TA effects also cause negative interactive SHAP values when TA is over the optimal temperature (∼25 • C).However, SWdif-VPD interactions cause GPP gains (positive interactive SHAP values) under high radiation and low VPD.High SWdif cause negative interactive SHAP values of SWdif-VPD at dry conditions due to plant stomatal closures.Among various PFTs, SWdif-TA and SWdif-VPD show similar pattern of responses for interactive SHAP values (figure S12).However, SWdif-SWC interactions present no significant responses among most PFTs except Shrub sites (figure S12).In Shrub sites, SWdif-SWC interactive SHAP values present positive responses to high SWC under more SWdif conditions, indicating higher diffuse fertilization efficiency with wetter soil conditions.These results are consistent with recent findings by Hemes et al (2020) though for wetland sites.With the longer timescales, our results of the interactions between SWdif and other climatic factors show no significant responses of interactive SHAP values at the daily or monthly timescales (figure S11), which may be related to limitations of sample numbers and smooth effects of climatic factors.
Previous studies found that SWC, VPD and CO 2 significantly influence terrestrial carbon uptake through water stress and CO 2 fertilization effects (Piao et al 2013, Yuan et al 2019).However, these two factors show much smaller impacts than temperature in the regulation of ecosystem carbon fluxes as revealed by the XGB-SHAP framework (figure 1).Such discrepancy is likely due to the application of allyear data in our analyses which magnifies the influences of TA through the phenological process.As a check, we used the same XGB-SHAP framework to explore the contributions of climatic factors to GPP variability during peak growing seasons across various timescales (figure S12).Without the regulation of phenology, we found the stronger impacts of SWC (12.4%),VPD (8.6%), and CO 2 (7.9%) than TA (6.2%) on GPP variability at the half-hourly timescale.However, even with such adjustment, we found the main feature remained the same with dominant impacts of SWdif on GPP variability at the short-term timescale.

Conclusions
Our results highlighted the important role of diffuse radiation in regulating the variability of ecosystem GPP.We found that SWdif dominates GPP variability of most PFTs (especially forest) at the hourly to sub-daily timescales.However, the contribution of SWdif decreases by about 8% from halfhourly to daily timescales, indicating that diurnal cycle of SWdif significantly affect the GPP variability.Even though temperature became more important at the daily or longer timescales, SWdif continues to act as a secondary driver in regulating GPP variability.Therefore, we suggest that the uncertainties of diffuse radiation should be constrained in the future projection of ecosystem carbon assimilation under various emission and warming scenarios.

Figure 1 .
Figure 1.SHAP values of climatic factors in the prediction of GPP at global FLUXNET sites.(a) Histogram of the mean absolute SHAP values of climatic factors at the half-hourly timescale.The errorbar indicates one standard deviation of global FLUXNET sites.(b) The bee swarm plots of GPP drivers.For each GPP driver, the range of y axis indicates sample size, and the color represents magnitude determined by the right colorbar.(c) Histogram of the mean absolute SHAP values from 3 hourly to monthly timescales.The GPP drivers include diffuse shortwave radiation (SWdif), air temperature (TA), direct shortwave radiation (SWdir), soil water content (SWC), vapor pressure deficit (VPD), carbon dioxide (CO2), surface pressure (P), wind speed (WS), and precipitation (PRE).

Figure 2 .
Figure 2. of mean absolute SHAP values for GPP variability caused by individual climatic factors for all and six plant functional types (PFTs) at (a) half-hourly and (b) daily timescales.The range of each boxplot indicates uncertainties and points indicate outliers for sites with the same PFT.Different colors represent varied climatic factors.Detailed site information is listed in tableS2.

Figure 3 .
Figure 3. Heatmaps of the mean absolute SHAP values for dominant GPP drivers and the interactive effects at (a) half-hourly and (b) daily timescales, and (c) their differences.The diagonal values indicate the main effects by single GPP drivers and non-diagnonal values indicate the interactive effects by the two drivers shown at x and y labels.Please note that the mean absolute SHAP values are calculated at global FLUXNET sites.

Figure 4 .
Figure 4.The main and interactive effects of climatic drivers on GPP for all and six PFTs across (a) half-hourly and (b) daily time scales, and (c) their differences.The diagonal values of the heatmap in (b) indicate the main effects by single GPP drivers, including SWdif, TA, SWdir, SWC and VPD.The non-diagonal values indicate the interactive effects by the two drivers shown at x and y axises.Each errorbar represents one standard deviation of the SHAP values at FLUXNET sites across all or individual PFTs.The color of histograms follows the heatmap in (b).

Figure 5 .
Figure 5. Summary of dominant climatic factors driving GPP variability for varied PFTs at different regions and timescales.The color of each square indicates one dominant climatic factor, and the number indicates one reference as listed in tableS3.
the noontime but show lower contributions at the sunrise and sunset periods (figureS4).The contributions by SWdir, SWC, VPD, and CO 2 show limited diurnal cycles with relatively high values in the noontime.For seasonal variations, both SWdif and TA show relatively higher contributions to GPP in summer and winter, and SWdir, SWC, VPD, and CO 2 make moderate contributions in summer (figure ations, SWdif shows high contributions at the sunrise, noontime, and sunset periods, indicating the stronger light controls of plant photosynthesis than other climatic factors.As the second GPP driver, TA becomes more important at