Spatial variations in terrestrial net ecosystem productivity and its local indicators

Multiple lines of evidence have demonstrated the persistence of global land carbon (C) sink during the past several decades. However, both annual net ecosystem productivity (NEP) and its inter-annual variation (IAV_NEP) keep varying over space. Thus, identifying local indicators for the spatially varying NEP and IAV_NEP is critical for locating the major and sustainable C sinks on land. Here, based on daily NEP observations from FLUXNET sites and large-scale estimates from an atmospheric-inversion product, we found a robust logarithmic correlation between annual NEP and seasonal carbon uptake–release ratio (i.e. U ∕ R). The cross-site variation in mean annual NEP could be logarithmically indicated by U ∕ R, while the spatial distribution of IAV_NEP was associated with the slope (i.e. β) of the logarithmic correlation between annual NEP and U ∕ R. Among biomes, for example, forests and croplands had the largest U ∕ R ratio (1.06 ± 0.83) and β (473 ± 112 g C m⁻² yr⁻¹), indicating the highest NEP and IAV_NEP in forests and croplands, respectively. We further showed that these two simple indicators could directly infer the spatial variations in NEP and IAV_NEP in global gridded NEP products. Overall, this study provides two simple local indicators for the intricate spatial variations in the strength and stability of land C sinks. These indicators could be helpful for locating the persistent terrestrial C sinks and provide valuable constraints for improving the simulation of land–atmospheric C exchanges.

1 Introduction

Terrestrial ecosystems reabsorb about one-quarter of anthropogenic CO₂ emission (Ciais et al., 2019) and are primarily responsible for the recent temporal fluctuations in the measured atmospheric-CO₂ growth rate (Randerson, 2013; Le Quéré et al., 2018). In addition, evidence based on eddy-flux measurements (Baldocchi et al., 2018; Rödenbeck et al., 2018), aircraft atmospheric budgets (Peylin et al., 2013) and process-based model simulations (Poulter et al., 2014; Ahlström et al., 2015) has shown a large spatial variability in net ecosystem productivity (NEP) on land. The elusive variation in terrestrial NEP over space refers to both of the substantial varying mean annual NEP and the divergent inter-annual variability (IAV) in NEP (i.e. IAV_NEP; usually quantified as the standard deviation of annual NEP) across space (Baldocchi et al., 2018; Marcolla et al., 2017). The mean annual NEP is related to the strength of carbon exchange of a specific ecosystem (Randerson et al., 2002; Luo and Weng, 2011; Jung et al., 2017), while IAV_NEP characterizes the stability of such a carbon exchange (Musavi et al., 2017). Thus, whether and how NEP and IAV_NEP change over space is important for predicting the future locations of carbon sinks on land (Yu et al., 2014; Niu et al., 2017).

Large spatial difference in terrestrial NEP has been reported from eddy-flux measurements, model outputs and atmospheric-inversion products. In addition, the global average IAV of NEP is large relative to global annual mean NEP (Baldocchi et al., 2018). More importantly, the spatial variations in NEP and IAV_NEP have been typically underestimated by the global flux-tower-based product and the process-based global models (Jung et al., 2020; Fu et al., 2019). These discrepancies have further revealed the necessity of identifying local indicators for the spatially varying NEP and IAV_NEP separately. The NEP in terrestrial ecosystems is determined by two components, including vegetation photosynthesis and ecosystem respiration (Reichstein et al., 2005), and their relative difference could determine the spatial variation in NEP (Baldocchi et al., 2015; Biederman et al., 2016). Many previous analyses have attributed the IAV_NEP at the site level to the different sensitivities of ecosystem photosynthesis and respiration to environmental drivers (Gilmanov et al., 2005; Reichstein et al., 2005) and biotic controls (Besnard et al., 2018; Musavi et al., 2017). For example, some studies have reported that IAV_NEP is more associated with variations in photosynthesis than carbon release (Ahlström et al., 2015; Novick et al., 2015; Li et al., 2017), whereas others have indicated that respiration is more sensitive to anomalous climate variability (Valentini et al., 2000; von Buttlar et al., 2018). However, despite previous efforts in a predictive understanding of the land–atmospheric C exchanges, the multi-model spread has not reduced over time (Arora et al., 2020). Therefore, it is imperative to explore the potential indicators for the spatially varying NEP, which could help attribute the spatial variation in NEP and IAV_NEP into different processes and provide valuable constraints for the global C cycle. Alternatively, the annual NEP of a given ecosystem can be also directly decomposed into net CO₂ uptake flux and CO₂ release flux (Gray et al., 2014), which are more direct components for NEP (Fu et al., 2019). It is still unclear how the ecosystem net CO₂ uptake and release fluxes would control the spatially varying NEP.

Conceptually, the total net CO₂ uptake flux (U) is determined by the length of the CO₂ uptake period (CUP) and the CO₂ uptake rate, while the total net CO₂ release flux (R) depends on the length of the CO₂ release period (CRP) and the CO₂ release rate (Fig. 1b). The variations in NEP thus could be attributed to these decomposed components. A strong spatial correlation between mean annual NEP and length of the CO₂ uptake period has been reported in evergreen needle- and broad-leaved forests (Churkina et al., 2005; Richardson et al., 2013; Keenan et al., 2014), whereas atmospheric-inversion data and vegetation photosynthesis models indicated a dominant role of the maximal carbon uptake rate (Fu et al., 2017; Zhou et al., 2017). However, the relative importance of these phenological and physiological indicators for the spatially varying NEP remains unclear.

Figure 1 Relationship between annual NEP and U ∕ R for 72 FLUXNET sites (of the form $\mathrm{NEP}=\mathit{\beta }\cdot \ln \left(\frac{U}{R}\right)$). (a) Dependence of annual NEP on the ratio between total CO₂ exchanges during net uptake (U) and release (R) periods (i.e. U ∕ R). Each line represents one flux site with at least 5 years of observations. (b) Conceptual figure for the decomposition framework introduced in this study. Annual NEP can be quantitatively decomposed into the following indicators: $\mathrm{NEP}=U-R$. (c) Distribution of the explanation of U ∕ R on temporal variability in NEP (R²) for FLUXNET sites.

In this study, we decomposed annual NEP into U and R and explored the local indicators for spatially varying NEP. Based on the eddy-covariance fluxes from the FLUXNET2015 dataset (Pastorello et al., 2017) and the atmospheric-inversion product (Rödenbeck et al., 2018), we examined the relationship between NEP and its direct components. In addition, we used the observations to evaluate the spatial variations in NEP and IAV_NEP in the FLUXCOM product and a process-based model (CLM4.5; Oleson et al., 2013). The major aim of this study is to explore whether there are useful local indicators for the spatially varying NEP and IAV_NEP in terrestrial ecosystems.

2 Materials and methods

2.1 Datasets

Daily NEP observations of eddy-covariance sites are obtained from the FLUXNET2015 Tier 1 dataset. The FLUXNET2015 dataset provides half-hourly data of carbon, water and energy fluxes at over 210 sites that are standardized and gap-filled (Pastorello et al., 2017). However, time series of most sites are still too short for the analysis of inter-annual variation in NEP. So only the sites that provided the availability of eddy-covariance flux measurements for at least 5 years are selected. This leads to a global dataset of 72 sites with different biomes across different climatic regions. Based on the biome classification from the International Geosphere–Biosphere Programme (IGBP) provided for the FLUXNET2015 sites, the selected sites include 35 forests (FORs), 15 grasslands (GRAs), 11 croplands (CROs), 4 wetlands (WETs), 2 shrublands (SHR) and 5 savannas (SAVs; Fig. S1 and Table S1 in the Supplement).

The Jena CarboScope inversion product combines high-precision measurements of atmospheric-CO₂ concentration with simulated atmospheric transport to infer the net CO₂ exchanges between land, ocean and atmosphere at large scales (Rödenbeck et al., 2018). Here, we used the daily land–atmosphere CO₂ fluxes from the s85_v4.1 version at a spatial resolution of $\mathrm{5}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$. Considering the relatively low spatial resolution of the Jena CarboScope inversion product, the daily fluxes were only used to calculate the local indicators for the spatially varying NEP at the global scale.

Daily NEP simulations from the Community Land Model version 4.5 (CLM4.5) were also used to calculate the local indicators for the spatially varying NEP at the corresponding flux tower sites. We ran the CLM4.5 model from 1985 to 2010 at a spatial resolution of 1^∘ with meteorological forcing from the Climate Research Unit and National Centers for Environmental Prediction (CRUNCEP). Here, NEP was derived as the difference between gross primary productivity (GPP) and total ecosystem respiration (TER), and TER was calculated as the sum of simulated autotrophic and heterotrophic respiration. The daily outputs from CLM4.5 were used to calculate the local indicators for the spatially varying NEP at both the global scale and the FLUXNET site level.

The FLUXCOM product presents an upscaling of carbon flux estimates from 224 flux tower sites based on multiple machine-learning algorithms and meteorological drivers (Jung et al., 2017). To be consistent with the meteorological forcing of the Jena CarboScope inversion product and the CLM4.5 model, we used the FLUXCOM CRUNCEPv6 products. In addition, in order to reduce the uncertainty caused by machine-learning methods, we averaged all the FLUXCOM CRUNCEPv6 products with different machine-learning methods. It should be noted that the inter-annual variability in the FLUXCOM product is driven by meteorological measurements and satellite data, which partially include information on vegetation state and other land surface properties. Daily outputs from FLUXCOM for the period 1985–2010 at 0.5^∘ spatial resolution were used to calculate the local indicators for the spatially varying NEP at both the global scale and the FLUXNET site level.

2.2 Decomposition of NEP and the calculations for its local indicators

The annual NEP of a given ecosystem can be defined numerically as the difference between the net CO₂ uptake and release (Fig. 1b). These components of NEP contain both photosynthesis and respiration flux, which directly indicate the net CO₂ exchange of an ecosystem. The total net CO₂ uptake flux (U) and the total net CO₂ release flux (R) can be further decomposed as

$$\begin{array}{}\text{(1)}& {\displaystyle }U=\stackrel{\mathrm{‾}}{U}\times \mathrm{CUP}\text{(2)}& {\displaystyle }R=\stackrel{\mathrm{‾}}{R}\times \mathrm{CRP},\end{array}$$

where CUP (d yr⁻¹) is the length of the CO₂ uptake period, and CRP (d yr⁻¹) is the length of the CO₂ release period; $\stackrel{\mathrm{‾}}{U}$ (g C m⁻² d⁻¹) is the mean daily net CO₂ uptake over CUP, and $\stackrel{\mathrm{‾}}{R}$ (g C m⁻² d⁻¹) represents the mean daily net CO₂ release over CRP. Many studies have reported that the vegetation net CO₂ uptake during the growing season and the non-growing season soil net CO₂ release are tightly correlated (Luo and Zhou, 2006; Zhao et al., 2016). Therefore, we further tested the relationship between annual NEP and U ∕ R (i.e. $\mathrm{NEP}\propto \frac{U}{R})$, which reflects the seasonal carbon uptake–release ratio. Consequently, NEP in any given ecosystem can be expressed as (Fig. S2)

$$\begin{array}{}\text{(3)}& \mathrm{NEP}=\mathit{\beta }\cdot \ln \left({\displaystyle \frac{U}{R}}\right),\end{array}$$

where the parameter β represents the slope of the linear relationship of $\mathrm{NEP}\propto \ln \left(\frac{U}{R}\right)$, indicating the site-level carbon uptake sensitivity. Based on the definitions of U and R, the ratio U ∕ R can be further written as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{U}{R}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}}\cdot {\displaystyle \frac{\mathrm{CUP}}{\mathrm{CRP}}}.\end{array}$$

The ratio of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ reflects the relative physiological difference between ecosystem CO₂ uptake and release strength, while the ratio of $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ is an indicator of net ecosystem CO₂ exchange phenology. Environmental changes may regulate these ecological processes and ultimately affect the ecosystem NEP. The slope β indicates the response sensitivity of NEP to the changes in phenology and physiological processes. All of β, $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ and $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ were then calculated from the selected eddy-covariance sites and the corresponding pixels of these sites in models. These derived indicators from eddy-covariance sites were then used to benchmark the results extracted from the same locations in models.

Figure 2 Relationship between annual NEP and U ∕ R for the Jena CarboScope inversion product (of the form $\mathrm{NEP}=\mathit{\beta }\cdot \ln \left(\frac{U}{R}\right)$). The black box indicates the location of the sample.

2.3 Calculation of the relative contributions

We further quantified the relative contributions of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ and $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ in driving the spatial variations in NEP:

$$\begin{array}{}\text{(5)}& \mathrm{NEP}=\mathit{\beta }\cdot \left[\ln \left({\displaystyle \frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}}\right)+\ln \left({\displaystyle \frac{\mathrm{CUP}}{\mathrm{CRP}}}\right)\right].\end{array}$$

For each eddy-covariance site, the parameter β was constant. Then, we used a relative-importance analysis method to quantify the relative contributions of these two ratios to the spatial variations in NEP. The algorithm was performed with the “relaimpo” package in R (R Development Core Team, 2011). The relaimpo package is based on variance decomposition for multiple-linear-regression models. We chose the most commonly used method named “Lindeman–Merenda–Gold” (LMG; Grömping, 2007) from the methods provided by the relaimpo package. This method allows us to quantify the contributions of explanatory variables in a multiple-linear-regression model. Across the 72 FLUXNET sites, we quantified the relative importance of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ and $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ to cross-site changes in NEP.

3 Results

3.1 The relationship between NEP and its direct components

To find local indicators for the spatially varying NEP in terrestrial ecosystems, we tested the relationship between NEP and its direct components (U and R) across the 72 flux tower sites. The results showed that annual NEP was closely related to the ratio U ∕ R (Fig. S2). The logarithmic correlations between annual NEP and U ∕ R were significant at all sites (Fig. 1a), and ∼90 % of R² values fall within a range of 0.7 to 1 (Fig. 1c).

In addition, the relationship between NEP and U ∕ R was also confirmed by the atmospheric-inversion product (i.e. Jena CarboScope inversion). The control of U ∕ R on annual NEP was robust in most global grid cells (i.e. $\mathrm{0.6}<R^{\mathrm{2}}<\mathrm{1}$). The coefficient of determination for this relationship was higher in 80 % of the regions but lower in North America (Fig. 2). These two datasets both showed that the indicator U ∕ R could successfully capture the variability in annual NEP.

3.2 Local indicators for spatially varying NEP

Across the 72 flux tower sites, the across-site variation in mean annual NEP was significantly correlated to mean annual $\ln \left(\frac{U}{R}\right)$ of each site (R²=0.65, P<0.01; Fig. 3a). In this network, the mean annual ratio $\ln \left(\frac{U}{R}\right)$ was a good indicator for cross-site variation in NEP. By contrast, the spatial variation in IAV_NEP was moderately explained by the slope (i.e. β) of the temporal correlation between NEP and $\ln \left(\frac{U}{R}\right)$ at each site (R²=0.39, P<0.01; Fig. 3b) rather than $\ln \left(\frac{U}{R}\right)$ (Fig. S3). The wide range of the ratio β reveals a large divergence in NEP sensitivity across biomes, ranging from 121 ± 118 g C m⁻² yr⁻¹ in shrubland to 473 ± 112 g C m⁻² yr⁻¹ in cropland.

Figure 3 Contributions of the two indicators in explaining the spatial patterns of mean annual NEP and IAV_NEP. (a) The relationship between annual mean NEP and $\ln \left(\frac{U}{R}\right)$ across FLUXNET sites (R²=0.65, P<0.01). The insets show the magnitude of $\ln \left(\frac{U}{R}\right)$ for different terrestrial biomes. (b) The explanation of β in IAV_NEP (R²=0.39, P<0.01). The insets show the distribution of parameter β for different terrestrial biomes. The number of years at each site is indicated with the size of the point.

The decomposition of indicator U ∕ R into $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ and $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ allowed us to quantify the relative importance of these two ratios in driving NEP variability. The linear regression and relative-importance analysis showed a more important role of $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ (58 %) than $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ (42 %) in explaining the cross-site variation in NEP (Fig. 4). Therefore, the spatial distribution of mean annual NEP was more strongly driven by the phenological changes.

Figure 4 The relative contributions of the local indicators in explaining the spatial patterns of mean annual NEP. (a) The linear regression between mean annual NEP with $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ (R²=0.33, P<0.01) and $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ (R²=0.25, P<0.01) across sites. (b) The relative contributions of each indicator to the spatial variation in mean annual NEP. The number of years at each site is indicated with the size of the point.

3.3 Simulated spatial variations in NEP by models

We further used these two simple indicators (i.e. $\frac{U}{R}$ and β) to evaluate the simulated spatial variations in NEP by the global flux-tower-based product (i.e. FLUXCOM) and a widely used process-based model at the FLUXNET site level (i.e. CLM4.5). We found that the low spatial variation in mean annual NEP in FLUXCOM and CLM4.5 could be inferred from their more converging $\ln \left(\frac{U}{R}\right)$ than flux tower measurements (Fig. 5). The underestimated variation in IAV_NEP in these modelling results was also clearly shown by the smaller β values (268.22, 126.00 and 145.08 for FLUXNET, FLUXCOM and CLM4.5, respectively; Fig. 5b).

Figure 5 Representations of the spatially varying NEP and its local indicators in the FLUXCOM product and the Community Land Model (CLM4.5) at the FLUXNET site level. (a) The variation in mean annual NEP and IAV_NEP is derived from FLUXNET, FLUXCOM and CLM4.5. Variation in mean annual NEP: the standard deviation of mean annual NEP across sites; variation in IAV_NEP: the standard deviation of IAV_NEP across sites. (b) Representations of the local indicators for NEP in FLUXNET, FLUXCOM and CLM4.5. The corresponding distributions of $\ln \left(\frac{U}{R}\right)$ and β are shown at the top and right. Significance of the relationship between annual NEP and $\ln \left(\frac{U}{R}\right)$ for each site is indicated by the circles. Closed circles: P<0.05; open circles: P>0.05. Note that the modelled results are from the pixels extracted from the same locations of the flux tower sites.

In addition, the spatial variations in NEP and IAV_NEP were associated with the spatial resolution of the product (Marcolla et al., 2017). Considering the scale mismatch between FLUXNET sites and the gridded product, we run the same analysis at the global scale based on the Jena CarboScope inversion product. At the global scale, the spatial variation in mean annual NEP can also be well indicated by $\ln \left(\frac{U}{R}\right)$ (Fig. 6). The larger net C uptake in FLUXCOM resulted from its higher simulations for $\ln \left(\frac{U}{R}\right)$. Furthermore, the larger spatial variation in IAV_NEP in CLM4.5 could be inferred from the indicator β.

Figure 6 Representations of the spatially varying NEP and its local indicators in the FLUXCOM product and the Community Land Model (CLM4.5) at the global scale. (a) The variation in mean annual NEP and IAV_NEP is derived from Jena CarboScope inversion, FLUXCOM and CLM4.5. Variation in mean annual NEP: the spatial variation in mean annual NEP; variation in IAV_NEP: the spatial variation in standard deviation in IAV_NEP. (b) Representations of the local indicators for NEP in Jena CarboScope inversion, FLUXCOM and CLM4.5.

4 Discussion

4.1 New perspective for locating the major and sustainable land C sinks

Large spatial differences in mean annual NEP and IAV_NEP have been well documented in previous studies (Jung et al., 2017; Marcolla et al., 2017; Fu et al., 2019). Here we provide a new perspective for quantifying the spatially varying NEP by tracing annual NEP into several local indicators. Therefore, these traceable indicators could provide useful constraints for predicting annual NEP, especially in areas without eddy-covariance towers.

Typically, the C sink capacity and its stability in a specific ecosystem are characterized separately (Keenan et al., 2014; Ahlström et al., 2015; Jung et al., 2017). Here we integrated NEP into two simple indicators that could directly locate the major and sustainable land C sink. Among biomes, forests and croplands had the largest $\ln \left(\frac{U}{R}\right)$ and β, indicating the strongest and the most unstable C sink in forests and croplands, respectively. However, the relatively lower β in shrublands and savannas should be interpreted cautiously. There are very few semi-arid ecosystems in the FLUXNET sites, while they represent a large portion of land at the global scale and have been shown to substantially control the inter-annual variability in NEP (Ahlström et al., 2015). The highest β implies that the land covered by cropland has the largest IAV_NEP. Therefore, the reported rapid global expansion of cropland may enlarge the fluctuations in land–atmosphere CO₂ exchange. In fact, the cropland expansion has been confirmed as one important driver of the recent increasing global vegetation growth peak (Huang et al., 2018) and atmospheric-CO₂ seasonal amplitude (Gray et al., 2014; Zeng et al., 2014).

4.2 Joint control of plant phenology and physiology on mean annual NEP

Recent studies have demonstrated that the spatio-temporal variations in terrestrial gross primary productivity are jointly controlled by plant phenology and physiology (Xia et al., 2015; Zhou et al., 2017). Here we demonstrated that the spatial difference in mean annual NEP was determined by both the phenology indicator $\frac{\mathrm{CUP}}{\mathrm{CRP}}$ (58 %) and the physiological indicator $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ (42 %). In addition, the lower contribution of the physiological indicator could partly be attributed to the convergence of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ across FLUXNET sites (Fig. S4).

The convergent $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ across sites was first discovered by Churkina et al. (2005) as 2.73 ± 1.08 across 28 sites, which included deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), and cropland and grassland. In this study, we found that the $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ across the 72 sites is 2.71 ± 1.61, which confirms the findings of Churkina et al. (2005). This spatial convergence of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ at the site level provides important constraints for global models that simulate large spatial variation in physiological processes (Peng et al., 2015; Xia et al., 2017). These findings imply that the phenology changes will greatly affect the locations of the terrestrial carbon sink by modifying the length of the carbon uptake period (Richardson et al., 2013; Keenan et al., 2014).

4.3 The simulated local indicators from gridded products

This study showed that the considerable spatial variations in mean annual NEP and IAV_NEP from global gridded products could also be inferred from their local indicators. The low variations in the U ∕ R ratio in CLM4.5 could be largely due to their simple representations of the diverse terrestrial plant communities in a few functional plant types with parameterized properties (Cui et al., 2019; Sakschewski et al., 2015). In addition, the higher U ∕ R ratio from the FLUXCOM product indicated its widely reported larger net C uptake (Fig. 6; Jung et al., 2020). Meanwhile, the ignorance of fire, land-use change and other disturbances could lead to the smaller β by allowing for only limited variations in phenological and physiological dynamics (Reichstein et al., 2014; Kunstler et al., 2016). Although the magnitude of IAV_NEP depends on the spatial resolution (Marcolla et al., 2017), we recommend that future model-benchmarking analyses use not only the global product compiled from the machine-learning method (Bonan et al., 2018) but also the site-level measurements or indicators (Xia et al., 2020).

4.4 Conclusions and further implications

In summary, this study highlights the changes in NEP and IAV_NEP over space on land and provides the U ∕ R ratio and β as two simple local indicators for their spatial variations. These indicators could be helpful for locating the persistent terrestrial C sinks where the $\ln \left(\frac{U}{R}\right)$ ratio is high but the β is low. Their estimates based on observations are also valuable for benchmarking and improving the simulation of land–atmospheric C exchanges in Earth system models. The findings in this study have some important implications for understanding the variation in NEP on land. First, forest ecosystems have the largest annual NEP due to the largest $\ln \left(\frac{U}{R}\right)$, while croplands show the highest IAV_NEP because of the highest β. Second, the spatial convergence of $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ suggests a tight linkage between plant growth and soil microbial activities in the non-growing season (Xia et al., 2014; Zhao et al., 2016). However, it remains unclear whether the inter-biome variation in $\frac{\stackrel{\mathrm{‾}}{U}}{\stackrel{\mathrm{‾}}{R}}$ is due to different plant–microbe interactions between biomes. Third, the within-site convergent but spatially varying β needs better understanding. Previous studies have shown that a rising standard deviation of ecosystem functions could indicate an impending ecological state transition (Carpenter and Brock, 2006; Scheffer et al., 2009). Thus, a sudden shift in the β value may be an important early-warning signal for the critical transition of carbon uptake sensitivity of an ecosystem. In this study, the atmospheric-inversion product shows low correlation between NEP and $\ln \left(\frac{U}{R}\right)$ in some boreal ecosystems, which might be due to the fact that the terrestrial NEP is not well constrained for these regions or because these boreal ecosystems are experiencing state transition. Therefore, the robustness in the relationship between annual NEP and $\ln \left(\frac{U}{R}\right)$ depends on the temporal stability of carbon uptake sensitivity for an ecosystem. In addition, the spatial variation in β reveals the differences in carbon uptake sensitivity across ecosystems. Furthermore, considering the limited eddy-covariance sites with long-term observations, these findings need further validation once the longer time series of measurements from more sites and vegetation types become available.