Climate and Vegetation Drivers of Terrestrial Carbon Fluxes: A Global Data Synthesis

The globally averaged air temperature has shown a warming trend since the industrial revolution (IPCC, 2013). Given that the terrestrial carbon (C) cycle plays an important role in global warming, it is critical to understand the key controls of the C metabolism of terrestrial ecosystems (Canadell et al., 2007). Over the past several decades, terrestrial ecosystems have generally been regarded as C sinks (Houghton, 2007), and the magnitude of terrestrial ecosystems in sequestering atmospheric carbon dioxide (CO₂) is of interest in mitigating climate changes because this terrestrial sink counters the anthropogenic increase in atmospheric CO₂ levels. (Friedlingstein et al., 2006) found that C-climate feedback was positive. Atmospheric CO₂ enters terrestrial ecosystems through photosynthesis, and terrestrial C is released back to the atmosphere through a variety of processes, collectively referred to as respiration (Trumbore, 2006).

Generally, the components of C flux include gross primary productivity (GPP), ecosystem respiration (ER), and net ecosystem production (NEP) (the difference between GPP and ER) (Reichstein et al., 2005). Eddy covariance is a method used in meteorology and other applications (e.g., micrometeorology and hydrology) to determine rates of gaseous exchange (Baldocchi, 2008). With the establishment of eddy covariance sites and relevant meteorological and biological measurements (Baldocchi et al., 2001; Valentini et al., 2000), a large body of data is available on the terrestrial ecosystem exchange of mass at stand level (Law et al., 2001; Wang et al., 2017). Measurements of GPP, ER and NEP are needed in order to determine the source-sink status of a specific ecosystem and to analyze how CO₂ exchange varies with the temporal variation in environmental conditions (Flanagan et al., 2002; Schmitt et al., 2010). Site-level observations have greatly advanced our understanding of patterns and environmental controls of terrestrial ecosystem C fluxes (Xu et al., 2014).

The temporal and spatial dynamics of GPP, ER and NEP are difficult to model or predict. Modeling efforts linking with observations at yearly scales are critical to progress in the field of terrestrial C fluxes. Analysis of global C flux data based on eddy covariance measurements makes it possible to identify ecosystem-scale CO₂ flux patterns that are not determined well in individual studies (Baldocchi, 2003; Baldocchi et al., 2015). The C flux measurements from individual studies have been synthesized to evaluate patterns of ecosystem C fluxes at regional (Valentini et al., 2000; Law et al., 2002; Kato and Tang, 2008; Chen et al., 2013; Yu et al., 2013; Xiao et al., 2013; Chen et al., 2019) and global scales (e.g., Luyssaert et al., 2007; Baldocchi, 2008; Yuan et al., 2009; Beer et al., 2010; Mahecha et al., 2010; Xu et al., 2014; Shao et al., 2015; Jung et al., 2017; Zhang et al., 2018). Although several studies have analyzed the spatial patterns of global C fluxes (Luyssaert et al., 2007; Migliavacca et al., 2011; Chen et al., 2015; Jung et al., 2017; Zhang et al., 2018), their datasets did not include all of the C flux data points worldwide and relevant vegetation characteristics [e.g., litterfall (LF), tree age (TA), diameter at breast height (DBH), tree height (TH), and basal area (BA)]. Despite hundreds of flux measurement sites, the availability of C flux data remains an important challenge for global synthesis. A large number of data points have been published in various literatures or are available publicly in the global FLUXNET (http://fluxnet.fluxdata.org/). Establishing such a dataset and synthesizing by using the data fusion technique can help in understanding the mechanisms of the C cycle and assessing the magnitude of the C budget (Wang et al., 2009).

It is likely that C fluxes vary across different biome types (e.g., non-wetland versus wetland), climate zones, and vegetation conditions. The influence of climate on C fluxes has been examined by several studies (e.g., Kato and Tang, 2008; Beer et al., 2010; Mahecha et al., 2010; Xiao et al., 2013; Chen et al., 2015). For instance, (Kato and Tang, 2008) indicated that CO₂ fluxes——particularly annual NEP——are determined mainly by mean annual temperature (MAT) and annual precipitation (AP). A global synthesis showed that GPP and ER are also related to MAT and AP (Chen et al., 2015). Our knowledge of ecosystem C cycling and its variability has been improved by investigating the global relationships between observed climate and site-level biosphere-atmosphere fluxes (Stoy et al., 2009). We have understood broadly how temperature and precipitation affect different C cycle processes. What we do have, however, is a quantitative understanding of the individual effects of temperature and precipitation and their comprehensive effects. Furthermore, the vegetation and environmental drivers of C fluxes are poorly understood and need to be further addressed. Particularly lacking, to the best of our knowledge, is compiled data focusing on synchronous measurements of NEP and vegetation characteristics and a comprehensive analysis of such data. Whether the dependency of C fluxes on climate differs with biome type, vegetation characteristics, or other factors is not well quantified. A better understanding of C flux dynamics will come from elucidating the integrated effects of climate and vegetation constraints on GPP, ER and NEP. Many more data on C cycling and vegetation characteristics in various biomes [e.g., forest, grassland (GL), wetland] make it possible to investigate the vegetation drivers of terrestrial C fluxes. In addition, the C cycle in wetlands is different from that in non-wetlands because of its special ecohydrology. The water layer in wetlands reduces soil CO₂ emissions, resulting in different C cycle patterns between wetlands and non-wetlands (Bridgham et al., 2006). Some studies that have synthesized flux observations from wetlands at regional scales typically used observations from a few sites (Turetsky et al., 2014; Petrescu et al., 2015). A recent article by (Lu et al., 2017) reported that ecosystem CO₂ fluxes of both inland and coastal wetlands are mainly regulated by MAT and AP, and the combined effects of MAT and AP can explain 71%, 54% and 57% of the variations in GPP, ER and NEP, respectively. However, the combined effects of MAT, AP, and vegetation characteristics [e.g., leaf area index (LAI)] on C fluxes have not been well explained, although many flux and vegetation measurements in wetlands worldwide have been well documented (e.g., Hirano et al., 2007; Hao et al., 2011; Beringer et al., 2013; Knox et al., 2015). The key questions that need to be addressed are whether the main drivers of C fluxes in wetlands are different from those in non-wetlands, and how to model C fluxes in wetlands using climate and vegetation characteristics.

In the present work, we assembled a comprehensive global dataset for terrestrial ecosystems that includes C fluxes (i.e., GPP, ER and NEP), biome categories, climate factors (temperature and precipitation), and vegetation traits (e.g., LAI). Although previous hypotheses have proposed that climate and vegetation characteristics may individually influence GPP, ER and NEP, the combined modeling of C fluxes on the basis of large-sample data on in-situ simultaneous measurements of climate and vegetation characteristics at the global scale remains lacking. Therefore, our present publicly available dataset is dedicated to re-examining previously hypothesized global relationships and developing new models. The objectives of the study were to: (2) present the dataset's structure and explain the global patterns of GPP, ER and NEP; (3) identify climate and vegetation drivers of these three fluxes for different biomes [e.g., forest, GL, cropland (CL), and wetland]; (3) determine the relative contributions of different drivers to the interannual and inter-site variabilities of these fluxes; and (4) establish models that fit the spatial and temporal variations in C fluxes in different biomes (wetland, forest, GL, and CL).

We compiled a dataset of 327 references from 296 sites that were useful for the synthesis analysis [Table S1 in electronic supplementary material A (ESM A)]. All of the references compiled in this study were published in peer-reviewed journals and covered C flux data from all continents except Antarctica. The oldest measurement year was 1991. The literature survey was intended to be as inclusive as possible up to 2016. Published studies that report annual NEP, GPP and ER were compiled using bibliographic databases including ISI Web of Science and Chinese Journal Net. The GPP, ER and NEP data compiled in this study were not directly downloaded from FLUXNET, because some of the eddy covariance data from FLUXNET have not been corrected. Instead, the C fluxes that appeared in the literature were used, but only where an appropriate correction method for the raw eddy covariance data was applied. We performed the data collection carefully to ensure the C fluxes (e.g., GPP, ER and NEP) and relevant environmental factors in the reviewed literature were used and that no part of them were removed. Studies were included if they met the following criteria: (2) they comprised annual NEP measurements only, i.e., studies of seasonal NEP measured during the growing season were not used; (3) the NEP field measurements were unmanipulated; and (3) flooded rice paddies, fenland, peatland, swamps, bogs, and marshes were considered as wetlands. Rice paddies were considered together with natural wetlands because the seasonal NEP patterns for flooded rice paddies are significantly different from upland fields, and land-use change from paddy rice cultivation to upland crop cultivation causes significant losses of C (Nishimura et al., 2008).

Measurement sites were distributed from 35^°39'S to 79^°56'N latitude and from 157^°25'W to 172^°45'E longitude, with most sites situated in the Northern Hemisphere (ESM A, Table S1). A total of 1194 measurements of annual NEP were collected from 296 sites, and most included annual GPP and ER measurement data. There are 870 data points for annual GPP and ER. The amount of site-year measurements included in the present dataset is more than that in previous studies (Luyssaert et al., 2007; Migliavacca et al., 2011; Chen et al., 2015). For example, the information on C fluxes, climate data and ecosystem properties (e.g., TA, TH and LAI) in the study of (Luyssaert et al., 2007) was only available for forest ecosystems. Our dataset includes six ecosystem types [CL, forest, GL, shrubland (SL), tundra (TD), and wetland]. In our dataset, most sites used the u* correction method to avoid underestimating the nighttime C effluxes (Lloyd and Taylor, 1994; Massman and Lee, 2002). Other methods, which was applied at relatively fewer sites, were used if the u* threshold was difficult to determine (e.g., Saitoh et al., 2005; Kato and Tang, 2008). We used the reported NEE values with the correction method that best considered the specific circumstances at each site.

For each study, whenever possible, we also noted the location (i.e., site name and country), latitude and longitude, measurement period, climate conditions, elevation, maximum C fluxes, soil and heterotrophic respiration (SR and HR), net primary productivity (NPP) and aboveground NPP (ANPP), vegetation characteristics, and soil properties (ESM A, Table S1). Variable categories such as climate conditions and vegetation characteristics are listed in Table 1. Vegetation characteristics included LF, fine root biomass (FR), NPP, ANPP, growing season length (GSL), TA, plant density of tree (PD), DBH, TH, BA, and maximum LAI. Soil properties included soil organic C storage (SOC), soil total nitrogen storage (STN), the ratio of C to N (C/N), pH, clay content (CLA), sand content (SAN), and bulk density (BD). Where MAT and AP were not available in the literature, they were derived from the Center for Climatic Research at the University of Delaware (http://climate.geog.udel.edu/~climate/html_pages/archive.html). The climate data were matched to the coordinates of the NEP study sites using a nearest-neighbor method. The maximum LAI (an important variable) data were compiled based only on studies that carried out manual measurements, rather than remote sensing, in order to avoid discrepancies between different datasets. Other potential variables (e.g., normalized difference vegetation index and greenness) were not compiled, as they are difficult to measure in-situ when C fluxes are measured at a specific site. The data were categorized into 10 biomes: broadleaf and needleleaf mixed forest (BNMF), CL, deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), GL, SL, tundra (TD), wetland [including forested wetland (FWL) and non-woody wetland (NWWL)] (Fig. 1). The FWL included mangroves and swamps, while the NWWL included marshland, mires, rice paddies, etc. The reason these 10 biome types were used was based on the consideration that leaf traits (i.e., broadleaf and needleleaf) may have a considerable influence on photosynthesis and respiration.

Figure 1.  The C flux sites compiled in this study.

We examined the distribution of annual C fluxes (i.e., GPP, ER and NEP) within and across biome types using the box-and-whisker plot. The box-and-whisker plot summarizes a data sample through five statistical measures: the minimum, lower quartile, median, upper quartile, and maximum. The mean values of each C flux were also computed for each biome type.

To clarify the issue as to how MAT and AP influenced the C flux component processes differently, for the growing season length (GSL) and the maximum GPP, ER and NEP, we grouped the global land ecosystems into 12 land climate classes. The scale for AP was divided into less than 0.4 m, 0.4-0.8 m, 0.8-1.5 m, and greater than 1.5 m, while that for MAT was less than 10^°C, 10^°C-20^°C, and greater than 20^°C. The AP of 0.4 m, 0.8 m and 1.5 m represents the thresholds between semiarid and semihumid zones, between semihumid and humid zones, and between humid and rainy zones, respectively. The 10^°C and 20^°C parts of the scale represent the temperature for active growth and vigorous growth of plants, respectively. The GPP, ER and NEP values were grouped into 12 land climate classes. The distribution of annual C fluxes (i.e., GPP, ER and NEP) in each climate class was determined using the box-and-whisker plot. Moreover, the GSL across different measurement sites and years could be explained as the MAT-related and AP-related functions, as the GSL was positively correlated with temperature in most climate zones and was closer to the seasonal precipitation pattern in tropical and subtropical regions (Liu et al., 2010; Ngeticha et al., 2014).

To identify the key factors controlling C fluxes, we conducted univariate regression analyses. In the univariate regression analyses, we examined the regulatory mechanisms of C fluxes by climate factors (MAT and AP), vegetation characteristics (e.g., LAI), and soil properties (e.g., SOC). The univariate regression analyses were used based on the following considerations. First, only potential factors influencing C fluxes with substantial determination coefficients (R²) and very low P values (P<0.001) in the univariate regression analyses can be considered as the variables in the following multivariate models. Second, the distribution patterns of C fluxes and potential drivers can be directly determined in the regression plots, which are necessary in choosing an appropriate mathematical expression in multivariate models. Third, MAT, AP and LAI have been commonly considered as independent predictive factors of C fluxes in previous univariate regression studies (Luyssaert et al., 2007; Lu et al., 2017; Hursh et al., 2017). Eventually, AP only explained ~10% (R²=0.108) of variations in LAI with a correlation coefficient (r) of 0.328 (ESM A, Table S1), while most of the variations in LAI (~90%) could not be explained by AP and may contribute to the variations in C fluxes. No significant correlation between MAT and LAI was found on the basis of available data (ESM A, Table S1).

Several aspects were considered in choosing the key drivers for the C flux models. First, the models for non-wetlands and wetlands might include different drivers. Second, both MAT and AP were considered as potential drivers in the model, as they have been found in previous investigations to be important controls of the variability in GPP and ER (Kato and Tang, 2008; Chen et al., 2013; Reichstein et al., 2013; Yu et al., 2013; Xiao et al., 2013). Third, only the potential vegetation and soil variables with enough data points and high significance can be included in the C flux models. When the mainly responsible drivers (climate and other site variables) were determined, they were all included in the model.

Before establishing the model, the relative contribution of each potential factor to each C flux was evaluated using the standardized coefficient in stepwise linear regression. To explain the dependence of annual C fluxes on MAT, AP, and one specific environmental control (X_i), we established the following model: \begin{eqnarray} \label{eq1} Y_{\rm c}&=&f_1({\rm MAT})f_2({\rm AP})f_3(X_i)\nonumber\\ &=&C_0e^{\alpha{\rm MAT}}\left(\frac{\rm AP}{{\rm AP}+\beta}\right)\left(\frac{X_i}{X_i+\gamma}\right). \ \ (1)\end{eqnarray} This model considers the interaction effects among the three variables of MAT, AP and LAI, whereas linear models cannot account for the interactions among variables. An exponential model has been commonly used to explain the relationship between C fluxes (GPP, ER and NEP) and temperature (Migliavacca et al., 2011; Ballantyne et al., 2017). A Michaelis-Menten function suggests declining positive effects of increasing precipitation or another environmental control (X_i) on C fluxes, which has been associated with the water and substrate supply on C assimilation and emission (Migliavacca et al., 2011; Exbrayat et al., 2013; Ballantyne et al., 2017). This model started from a widely used climate-driven model proposed by (Raich et al., 2002) and further modified by (Reichstein et al., 2003) in the form of a climate-and-LAI-driven model. In the equation, Y_c denotes a specific C flux item (e.g., GPP); C₀ is the Y_c at 0^°C without precipitation and X_i (one specific environmental control) limitations; α (^°C^-1) determines the increasing rate of Y_c with the increase of MAT; and β (m) and γ represent the half-saturation constants of the Michaelis-Menten function determining the relationship between Y_c and precipitation or X_i, respectively. The relationship between precipitation (or X_i) and observed Y _c can be approximated using a Michaelis-Menten function, with increasing rates of precipitation (or X_i) having lesser impacts on fluxes. The model's parameters (C₀, α, β and γ) were estimated using a least sum of residual squares method. The modeling errors were determined by performing a bootstrapping procedure using resampling from the original dataset to create 100 different datasets (Cameron et al., 2008; Keenan et al., 2013). The GPP, ER and NEP within each of the 12 climate classes were modeled by the climate (MAT and AP) and vegetation variable (LAI) within each climate class. The model structure was the same as in Eq. (2).

Six statistical measures, including the adjusted coefficient of determination (R²), probability of obtaining a test statistic (P), root-mean-square error (RMSE), model efficiency (ME), mean absolute error (MAE) (Janssen and Heuberger, 1995), and Akaike information criterion (AIC) (Akaike, 1974), were used to evaluate the performance of the established models. The RMSE, ME, MAE and AIC were computed by the following equations: \begin{eqnarray} \label{eq2} {\rm RMSE}&\!\!=\!\!&\sqrt {\frac{\sum_{i=1}^{n}(Y_{{\rm c},{\rm MOD},i}-Y_{{\rm c},{\rm OBS},i})^{2}}{n}} ; \ \ (2) \\ \label{eq3} {\rm ME}&\!\!=\!\!&\frac{[\sum_{i=1}^{n}\!(Y_{{\rm c},{\rm OBS},i}\!\!-\!\!\overline{Y_{{\rm c},{\rm OBS}}})^{2}\!-\!\sum_{i=1}^{n}\!(Y_{{\rm c},{\rm MOD},i}\!-\!Y_{{\rm c},{\rm OBS},i})^{2}]} {[\sum_{i=1}^{n}(Y_{{\rm c},{\rm OBS},i}\!-\!\overline{Y_{{\rm c},{\rm OBS}}})^{2}]}\,; \ \ (3) \\ {\rm MAE}&\!\!=\!\!&\frac{\sum_{i=1}^n|Y_{{\rm c},{\rm MOD},i}-Y_{{\rm c},{\rm OBS},i}|}{n} ; \ \ (4) \\ \label{eq4} {\rm AIC}&\!\!=\!\!&2k-2\ln (L) . \ \ (5)\end{eqnarray}

In above equations, Y _{c, MOD} and Y _{c, OBS} are the modeled and observed values of C fluxes, respectively; $\overline {Y_\rm c,\rm OBS}$ denotes the mean of Y _{c, OBS,i}; n is the number of samples; k defines the number of parameters in the statistical model; and L is the maximized value of the likelihood function for the estimated model.

All statistical analyses were performed with SPSS19.0 for Windows. All levels of significance reported are P<0.05 unless otherwise noted.

The magnitude of annual C fluxes varied across biomes (Figs. 2a-c). Annual GPP and ER of terrestrial ecosystems ranged from 3.417 to 305.500 mol C m^-2 yr^-1 and from 2.833 to 325.583 mol C m^-2 yr^-1, respectively. On average, EBF had the highest annual GPP and ER, while TD had the lowest annual GPP and ER. The highest annual GPP and ER were from Simpang Pertang (Malaysia) and Palangka Raya (Indonesia), respectively, while the lowest annual GPP and ER were from Inner Mongolia (China) and Daring Lake (Canada), respectively. The range of annual GPP and ER was greatest for wetlands, as they included peatland in the tropical zone and fenland in the frigid zone (ESMA, Table S1). For each biome, the mean annual GPP was slightly higher than the mean annual ER. Annual NEP varied from -113.333 to 116.667 mol C m^-2 yr^-1 (Fig. 2c). Compared with GPP and ER, annual NEP showed lower variance across biomes. Mean annual NEP was 20.417, 17.167, 22.417, 19.917, 25.833, 18.750, 8.083, 4.833, 2.833, 12.538 and 7.065 mol C m^-2 yr^-1 for the BNMF, CL, BDF, DNF, EBF, ENF, GL, SL, TD, FWL and NWWL biomes, respectively. Among the biomes, the highest and lowest mean annual NEP was for EBF and TD, respectively. More than 80% NEP values were greater than zero, indicating most examined ecosystems across sites and years were C sinks. The present GPP and ER dataset for different forest biomes also validated the theory that broadleaf forest has higher photosynthesis and respiration rates than needleleaf forest. The net C assimilation (NEP) was also slightly higher for broadleaf forest than needleleaf forest.

Figure 2.  Box-and-whisker plots for C fluxes: (a) global GPP; (b) ER; (c) NEP.

As indicated in the box-and-whisker plots (Fig. 3), the distribution of annual NEP showed clear patterns with the changing climate classes. These results indicated that MAT and MAP influence the C flux component processes differently. The median of NEP declined with the increase in the MAT scale when AP was scarce (<0.4 m). However, NEP increased with the increase in MAT when AP ranged from 0.4 to 1.5 m. The maximum NEP appeared in the MAT range of 10^°C-20^°C when AP was plentiful (>1.5 m). Across the 12 climate classes, the maximum NEP appeared at the upper end of the MAT scale (MAT>20^°C) and in the moderate part of the AP scale (0.8 m < AP < 1.5 m). These patterns of NEP depended on the variations in GPP and ER. Both GPP and ER increased with the increase in MAT and AP, with the maximum GPP and ER at the upper end of the MAT scale (MAT>20^°C) and AP scale (AP>1.5 m). For all levels of precipitation zones (<0.4 m; 0.4-0.8 m; 0.8-1.5 m; >1.5 m), the median of GPP generally increased with the increase in MAT [Table S1 in electronic supplementary material B (ESM B)]. A similar phenomenon could be seen for the median of ER (ESM B, Table S1), although the tendency was not as obvious as for GPP.

Figure 3.  Box-and-whisker plots for C fluxes at different MAT and AP scales. The median values of NEP are shown above the NEP plot.

Annual GPP exhibited strong relationships with climate factors. For all biomes, GPP was positively related to MAT and AP (Figs. 4a and b). About 17.6% and 32.7% of the variations in GPP could be explained by MAT and AP, respectively. The relationship between GPP and MAT was fitted by an exponential equation. GPP increased exponentially with the increase in MAT for both non-wetland and wetland (i.e., FWL and NWWL) biomes (Figs. 4d and g). Compared with temperature, precipitation had a stronger correlation with GPP for non-wetland biomes. AP explained 42.9% of the variability in GPP for non-wetland biomes (Fig. 4e).

Figure 4.  Relationship between GPP and environmental variables: (a-c) MAT, AP and LAI, respectively, for all biomes; (d-f) MAT, AP and LAI, respectively, for non-wetland sites; (g-i) MAT, AP and LAI, respectively, for wetland sites. The P values in all panels are less than 0.001.

GPP was also strongly correlated with LAI (Fig. 4c). LAI was an important factor determining GPP for both non-wetland and wetland (i.e., FWL and NWWL) biomes, explaining 34.5% (P<0.001) and 74.6% (P<0.001) of the variability, respectively. For all biomes, LAI explained 40.2% of the variability in GPP. For wetland (i.e., FWL and NWWL), GPP was significantly correlated to MAT and LAI (Figs. 4g and i); GPP was not significantly correlated with AP in this ecosystem (Fig. 4h).

Similar to GPP, ER responded to a suite of drivers, including temperature, precipitation and vegetation productivity (Fig. 5). In the univariate regressions, MAT, AP and LAI were the site-specific variables that correlated most strongly with ER. For all biomes, 19.1%, 32.4% and 33.2% of the variations in ER could be explained by MAT, AP and LAI, respectively (Figs. 5a-c). Ecosystems with higher temperature, precipitation and LAI had a larger value of ER as C output. For wetland (i.e., FWL and NWWL), the drivers of ER were MAT and LAI (Figs. 5g and i); ER did not exhibit clear precipitation dependency in this ecosystem (Fig 5h). Similar to GPP, ER was also not significantly correlated with soil properties (ESM A, Table S1). The relation between annual ER and GPP is shown in Fig. 6. The slope of the relation was 0.805 for all biomes, 0.767 for non-wetland, and 1.140 for wetland (i.e., FWL and NWWL) (intercepts of 4.903, 7.847 and -11.498 mol C m^-2 per year, respectively). The ER/GPP ratio averaged 0.805 across all sites and years, with values greater than 1 (i.e., net loss of CO₂) for wetland (i.e., FWL and NWWL).

Figure 5.  Relationship between ER and environmental variables: (a-c) MAT, AP and LAI, respectively, in all biomes; (d-f) MAT, AP and LAI, respectively, for non-wetland sites; (g-i) MAT, AP and LAI, respectively, for wetland sites. The P values in all panels are less than 0.001.

Figure 6.  Relationship between ER and GPP: (a) all biomes; (b) non-wetland; (c) wetland. The P values in all panels are less than 0.001.

Annual NEP was weakly correlated with climate factors (i.e., MAT and AP) (Figs. 7a, b, d, e, g and h). Similar to GPP and ER, NEP was significantly (P<0.001) correlated with LAI when all biomes (Fig. 7c) and non-wetland (Fig. 7f) were examined. No significant correlation between NEP and LAI was found for wetland (i.e., FWL and NWWL) (Fig. 7i). Although the R² for the function explaining the relationship between NEP and LAI was relatively low (R²=0.106), the number of independent data points for the fitting was large, ensuring the reliability of the regression function. Moreover, annual NEP was significantly correlated with GPP (ESM B, Fig. S1). The slope of the relation was 0.195 for all biomes, 0.234 for non-wetland, and -0.140 for wetland (i.e., FWL and NWWL). Annual NEP increased with the increase in GPP for non-wetland, and decreased with the increase in GPP for wetland (i.e., FWL and NWWL).

Figure 7.  Relationship between NEP and environmental variables: (a-c) MAT, AP and LAI, respectively, in all biomes; (d-f) MAT, AP and LAI, respectively, for non-wetland sites; (g-i) MAT, AP and LAI, respectively, for wetland sites. The P values in (a-f) are less than 0.001; the P values in (g-i) are 0.003, 0.019, and 0.921, respectively.

NEP in our dataset varied with measurement period across the global 24-year C fluxes measurements (ESM B, Fig. S2). A quadratic function could explain the relationship between GPP (or NEP) and measurement period, suggesting the interannual variability of the GPP (or NEP). NEP during the 1993-97 and 2009-14 time periods was less than that during other time periods (ESM B, Fig. S3).

NEP at some multiple-year flux sites showed considerable variability when the interannual variation of NEP was considered. NEP could be quite high in the first few years after a long dry period, and then became quite low or neutral even though MAT and AP did not change much for a few years after the long dry period (ESM A, Table S1). The representative sites that had such phenomenon were Qianyanzhou (China), Duke forest (USA), Turkey Point (Canada), Pellston (USA), Prince Albert (Canada), Hyytiälä (Finland), Albuquerque (USA), Dresden (Germany), and Twitchell (USA) (ESM A, Table S1). Plants might show marked responses to increasing precipitation after a long dry period, which could be explained as priming effects. After the dry and humid period, the priming effects of plants to changing precipitation might not be obvious, because plants have acclimatized to the precipitation fluctuation.

The relative contribution of MAT, AP and LAI to each C flux is indicated in Table S2 (ESM B). For the non-wetland, the relative contribution of AP to GPP/ER/NEP in the climate factor-based model was higher than that of MAT, while the relative contribution of LAI to GPP/ER/NEP in the climate- and LAI-based model was higher than that of both MAT and AP. For the wetland (i.e., FWL and NWWL), the relative contribution of MAT to GPP/ER/NEP was higher than that of both AP and LAI.

The models including climate factors and LAI are shown in Table 2. For all biomes, the climate factor model named MAT&AP-model explained 38.3% of variations in GPP. The combined contribution of MAT and AP to the spatial and temporal variations in GPP increased significantly compared with the single-factor contribution of MAT or AP (Fig. 4 and Table 2). For all biomes, the RMSE, ME, and MAE for the climate-based GPP model were 46.265 mol C m^-2 yr^-1, 0.383, and 35.993 mol C m^-2 yr^-1, respectively. For non-wetland, the climate factor model explained 42.9% of variations in GPP. Statistical analyses indicated that the climate- and LAI-based model named MAT&AP&LAI-model performed well in simulating the temporal and spatial variability in GPP across various biomes (Table 2). The model fit of GPP of the climate-based equation related to MAT and AP was improved when the LAI was incorporated into the model, suggesting that GPP increases with increasing LAI. For all biomes, MAT&AP&LAI-model appeared to be the best fit for the variations in GPP when the statistical measures of n, R², P, RMSE, ME, MAE and AIC were comprehensively evaluated (Table 2). Moreover, for FWL and NWWL, the fit of MAT&LAI-model for GPP had high values of R² (0.914) and ME (0.887) and low values of RMSE (23.387) and MAE (16.688), suggesting that they strongly explain the variations in GPP.

The climate factor model named MAT&AP-model explained 36.4% of variations in ER. The RMSE, ME, and MAE for the climate-based ER model were 46.265 mol C m^-2 yr^-1, 0.364, and 30.904 mol C m^-2 yr^-1, respectively. For non-wetland, the climate factor model explained 45.4% of variations in ER. By integrating the interaction between climate and LAI, the explanation of the regression equation with respect to ER also increased in comparison with MAT&AP-model, when the statistical measures of n, R², P, RMSE, ME, MAE and AIC were comprehensively evaluated (Table 2). Moreover, for FWL and NWWL, MAT&AP-model explained 92.9% (R²=0.929) of the variations in ER, with low values of RMSE (21.151 mol C m^-2 yr^-1) and MAE (14.811 mol C m^-2 yr^-1).

Compared with GPP and ER, MAT&AP-model explained a lower proportion of the variability in NEP for FWL and NWWL. For all biomes, the RMSE, ME, and MAE for the climate-based NEP model were 23.042 mol C m^-2 yr^-1, 0.156, and 16.467 mol C m^-2 yr^-1, respectively (Table 2). The climate factor model explained 7.7% and 14.6% of the variations in NEP for non-wetland and wetland (i.e., FWL and NWWL).

To explore the effects of land use or disturbance (i.e., TA, TH, DBH and BA) on NEP, we further examined the relationship between NEP and these variables (Figs. 8a-d). As shown in Fig. 8a, NEP increased significantly (P<0.001) with the increase in TA when TA was less than 15 yr. However, NEP declined significantly (P<0.001) with the increase in TA when TA was more than 15 yr, indicating that young trees aged from 10 to 20 yr can absorb more C from the atmosphere. Similar to TA, the trees with height ranging from 10 to 20 m had the highest ability to absorb C. The relationship between NEP and DBH could be well explained by a quadratic model, with a determination coefficient of 0.185.

Figure 8.  Relationship between NEP and tree characteristics (TA, TH, DBH and BA) due to land use or disturbance: (a) TA; (b) TH; (c) DBH; (d) BA. The several black outliers are not used for establishing the regression line in panel (d). The P values in all panels are less than 0.001.

Moreover, with over 1000 annual observations for each of those three fluxes from about 300 stations, a number of sites had multiple-year observations. In our analysis, multiple observations from one station were treated as independent observations. Therefore, both inter-site and interannual variations in C fluxes could be determined and modeled. If the multiple-year observations from one station were considered as one sample and modeled on the basis of the averaged multiple-year fluxes and multiple-year environmental variables at each site, the relationship between C fluxes and environmental variables on the basis of integrated site-year data (ESM B, Fig. S4) was similar to that on the basis of all inter-site and interannual data (Figs. 4, 5 and 7).

AP was the most important factor controlling annual GPP and ER across non-wetland biomes. Compared with temperature, precipitation explained a higher proportion of the variance in GPP and ER. The univariate regressions indicated that precipitation explained 42.9% of the variations in GPP and 44.6% of the variations in ER for non-wetland (Figs. 4e and 5e). Temperature was another important climatic factor controlling the variability in GPP and ER. Overall, warmer and wetter sites had higher GPP and ER, while drier and/or colder sites had lower GPP and ER.

The relationship between GPP (or ER) and climate factors has been reported by a number of regional and global studies (Law et al., 2002; Kato and Tang, 2008; Chen et al., 2013; Reichstein et al., 2013; Yu et al., 2013; Xiao et al., 2013; Xu et al., 2014). For instance, (Chen et al., 2013) reported that the combined effects of MAT and AP accounted for 85% and 81% of the spatial variations in GPP and ER in the Asian region, respectively. (Xu et al., 2014) examined the global patterns and environmental controls of forest C balance, and found that both production and respiration increased with MAT and exhibited unimodal patterns along a gradient of precipitation. For wetland (i.e., FWL and NWWL), only MAT was responsible for the spatial and temporal variations in GPP and ER. AP was not significantly correlated with annual GPP and ER, as water is not the common limitation to vegetation growth in wetland.

Our study showed that NEP was significantly correlated with MAT and AP (Fig. 7), which was in agreement with the results of some regional studies showing that higher air temperature and/or precipitation can lead to higher net C uptake (Yu et al., 2013; Xiao et al., 2013). (Chen et al., 2013) reported that MAT and AP explained 36% of the variations in NEP in Asian ecosystems. Annual NEP was also found to be positively related to annual temperature and precipitation in Chinese terrestrial ecosystems (Xiao et al., 2013). Moreover, studies have shown that climate variation can be directly responsible for short- but not long-term variation in forest-atmosphere C exchange (Richardson et al., 2007). At global scales, our study showed that explaining NEP variability only via climatic variables is ineffective.

In this study, precipitation was used as a proxy of water availability. Although precipitation is relevant, other proxies (e.g., evapotranspiration, dryness, soil moisture) are also relevant. (Yi et al., 2010) showed that NEP was generally a function of temperature at colder sites, relative to a function of dryness at warmer sites. This suggests that NEP is controlled by a complex interaction of climate factors (i.e., temperature and dryness). Models incorporating soil moisture perform differently to models using precipitation when simulating soil C fluxes (McGuire et al., 2000; Exbrayat et al., 2013; Chen et al., 2014; Hursh et al., 2017). Although factors relevant to the relationship between precipitation, dryness and soil moisture are difficult to disentangle at broad scales, temperature and soil texture are thought to be the potential key drivers (Davidson et al., 2000). Addressing this relationship requires more measurement data across a larger number of site-years. Long-term observations of GPP, ER and NEP in various climate zones are needed to verify their variability at different time scales and in different regions.

Vegetation drivers have been shown to be critical to the inter-site and interannual variations in C fluxes (Jung et al., 2011). In the present study, LAI was an important factor determining GPP and ER for both non-wetland and wetland (i.e., FWL and NWWL) biomes. Previously, few investigations have reported the correlation between annual GPP and LAI (Luyssaert et al., 2007). On an annual basis, (Law et al., 2002) found a poor correlation between GPP and LAI, possibly because the amount of data points was limited in their study. Indeed, LAI is a unique biophysical factor accounting for the differences in phenological development, assimilation and biomass growth in plant canopies (Schmitt et al., 2010). Leaf area exerts a major influence on canopy photosynthesis (Tappeiner and Cernusca, 1996; Saigusa et al., 1998; Restrepo-Coupe et al., 2013; Wu et al., 2016; Baldocchi et al., 2018), which also provides assimilates for the respiration of roots and soil microorganisms (Reichstein et al., 2003; Bahn et al., 2008, Bahn et al., 2009; Bond-Lamberty and Thomson, 2010). By testing the pairwise relationship between ER and different site characteristics, (Migliavacca et al., 2011) found that the ecosystem LAI showed the closest correlation with ER, and thus LAI was the best explanatory variable of the ER variability.

Besides LAI, modeling and in-situ measurement studies have highlighted the importance of biotic drivers, such as stand age, to spatial patterns of C fluxes (Migliavacca et al., 2011). Substrate supply is thought to be more important than temperature in determining ER across European forests (Janssens et al., 2001). The endogenous processes regulate a large part of the daily variation in terrestrial NEP (de Dios et al., 2012). Moreover, ecosystems under the climate optimum often reach their photosynthetic and respiratory potentials, resulting in the larger contribution of biological effects (compared with climate effects) to the interannual variability in NEP (Shao et al., 2015). A number of previous studies have shown that NEP is affected by interactions of microclimate, canopy structure, and the photosynthetic and respiratory physiology of plants (Flanagan and Johnson, 2005; Luyssaert et al., 2007; Schmitt et al., 2010), which is essentially in agreement with our results.

MAT&AP-model, MAT&AP&LAI-model, and MAT& LAI-model relied on statistical relations between annual GP/ER/NEP and some key controls (climate or LAI), rather than a complete understanding of the mechanisms involved. It was necessary that these simulations generated the correct order of magnitude of measured values and the general patterns from low to high C fluxes. Although previous studies have employed climate variables to simulate C fluxes (Schmitt et al., 2010; Migliavacca et al., 2011; Yu et al., 2013; Chen et al., 2015), our models included the most inclusive data points of C fluxes and relevant climate and vegetation variables, increasing the simulation efficiency. These results obtained with the bootstrap estimation of MAT&AP-model and MAT&AP&LAI-Model indicated that the developed model described the GPP and ER for all biomes well.

The derived parameterization of MAT&AP&LAI-model recorded in Table 2 may be considered as an optimized parameterization for the application of the model at global scales. One of the main advances introduced by this model formulation is the incorporation of LAI as a driver of the GPP and ER. This variable was necessary to improve the description of both the temporal and spatial dynamics of GPP and ER. In this study, wetland (i.e., FWL and NWWL) was insensitive to precipitation, possibility owing to the presence of water. For wetland, MAT&LAI-model explained more than 90% of variations in GPP and ER, suggesting MAT and LAI are the two main drivers of GPP and ER.

However, the models shown in Table 2 are relatively crude and did not vary by more specific biome types (i.e., between EBF, ENF, GL, TD, etc.), in spite of the significant variation in biotic controls that would be expected. The unexplained variability in the multiple regressions indicated that GPP, ER and NEP are probably related to a number of factors, including species composition, vegetation characteristics, and spatial variability in nutrient availability (Bahn et al., 2006; Schmitt et al., 2010). The information about vegetation characteristics other than LAI was limited, making it impossible to incorporate these factors into the model. Moreover, the lack of vegetation data for each site and each year increased the difficulties of biome-specific modeling. More sophisticated analyses are expected if more such vegetation data become available.

The NEP was relatively constant across all types of biomes (Fig. 2c). The low variability in NEP constituted a major reason why its measurement and modeling remained so difficult. By contrast, GPP and ER showed great variability across biomes, and even across sites and years for a specific biome (Figs. 2a and b); the variability in GPP and ER was generally in accordance with that in temperature and precipitation. NEP responded to climate drivers with relatively low R² (Fig. 7). Each biome (from EBF in the tropical area to TD in the frigid zone) can act as a net C sink (i.e., NEP>0), while controls on the balance between photosynthesis and respiration varied with sites and time scales. The weather conditions like drought, the phenology of plants, and the supply of decomposable litter or exudates may create temporary C storage or loss in ecosystems (Trumbore, 2006). However, we did find the climate and vegetation controls on NEP (Figs. 7 and 8), providing the clues for modeling NEP. The simulation efficiency of NEP may be improved if available site-specific vegetation data are enough for a particular biome in the future.

The models simulating the variations in the C fluxes of GPP, ER and NEP varied at the different scales of MAT and AP (ESM B, Table S3). LAI was a key driver of C fluxes of GPP, ER and NEP when the MAT was at the scale of less than 10^°C. The influence of MAT on C fluxes was significant at the scale of lowest MAT (MAT<10^°C) and least AP (AP<0.4 m). AP was one of the key drivers of C fluxes of GPP, ER and NEP from the semiarid to -rainy zones (0.4 m< AP<1.5 m) when the MAT was above 20^°C. The LAI was the only driver of GPP and ER when the MAT and AP was greater than 20^°C and 1.5 m, respectively. Overall, the variations in GPP, ER and NEP were mainly controlled by LAI (a key indicator of plant growth) when the MAT and AP were relatively low or high, while their variations were controlled by MAT, AP or LAI at the moderate scale of MAT and AP. The distribution of C fluxes and modeling in 12 land climate classes would not only help to determine how temperature, precipitation and vegetation influence C flux components differently, but also provide new clues as to how to improve world climate classification (e.g., Köppen-Geiger climate classification) (Rubel and Kottek, 2010) using C flux patterns (Fig. 3). A new climate classification scheme in the future on the basis of C flux patterns may be more effective as it considers the plant growth and C assimilation. Further subdivision of the climate classes (e.g., dividing the semiarid zone into AP less than 0.2 m and within the range of 0.2 to 0.4 m) is also needed by compiling more field measurement data.

Further analysis using Pearson's correlation coefficient suggested that maximum C fluxes of GPP_max, ER_max and NEP_max over the growing season were significantly (P<0.001) correlated with climate factors (MAT and AP) (ESM B, Table S4). The variations in GPP_max and NEP_max were controlled by LAI, while the variations in ER_max were not significantly (P>0.05) correlated with LAI. MAT, AP and LAI significantly (P<0.001) influenced the variations in GSL, and GSL was significantly (P<0.001) correlated with annual C fluxes (GPP, ER and NEP). Overall, the maximum C fluxes, annual C fluxes, GSL and climate factors had close relationships.

In our study, ER was strongly correlated with GPP, which could be validated by a number of studies at different spatial and temporal scales (Lasslop et al., 2010). (Larsen et al., 2007) observed that ER depended strongly on photosynthesis in a temperate heath. (Janssens et al., 2001) and (Reichstein et al., 2007) reported that the annual ER increased linearly with GPP across European forests. (Yu et al., 2013) found that the spatial patterns of GPP and ER of all terrestrial ecosystems in China showed a positive correlation. The spatial and temporal variations in ER and GPP were also coupled when this analysis was expanded to global ecosystems (Chen et al., 2015). Overall, this tight coupling confirms that the variation in photosynthate availability is the dominant driver of respiration and the coupling of GPP and ER should be taken into account in terrestrial C cycle models.

Our study provides new insights into the drivers of the interannual and inter-site variability of NEP, GPP and ER. Across terrestrial biomes, climate factors were found to play an important role in determining the GPP and ER. For non-wetland, GPP and ER tended to increase towards a warmer and wetter environment. For wetland (i.e., FWL and NWWL), higher GPP and ER appeared in warmer conditions, while precipitation was not the factor limiting site productivity and respiration. Meanwhile, the climate roles were mediated by the vegetation characteristics, such as LAI. However, climate only showed little impact on the variability in NEP, which was relatively constant across all types of biomes. In addition, the vegetation controls (e.g., LAI and TA) on NEP provided the clues to model NEP. Overall, the effects of LAI on GPP, ER and NEP indicated that canopy-level measurement is critical for accurately estimating ecosystem-atmosphere exchange of CO₂. Climate and LAI can influence C flux components differently in different climate regions. This synthesis study highlights that the responses of ecosystem-atmosphere exchange of CO₂ to climate and vegetation variations are complex, which poses great challenges to models seeking to represent terrestrial ecosystem responses to climatic variation.

Electronic supplementary material. Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-019-8194-y.