# Regional and seasonal partitioning of water and temperature controls on global land carbon uptake variability

Kai Wang, Ana Bastos, Philippe Ciais, Xuhui Wang, Christian Rödenbeck, Pierre Gentine, Frédéric Chevallier, Vincent W. Humphrey, Chris Huntingford, Michael O'Sullivan, et al.

► To cite this version:

Kai Wang, Ana Bastos, Philippe Ciais, Xuhui Wang, Christian Rödenbeck, et al. Regional and seasonal partitioning of water and temperature controls on global land carbon uptake variability. Nature Communications, 2022, 13, 10.1038/s41467-022-31175-w. insu-03824370

HAL Id: insu-03824370

https://insu.hal.science/insu-03824370

Submitted on 21 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-022-31175-w

<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# Regional and seasonal partitioning of water and temperature controls on global land carbon uptake variability

Kai Wang<sup>1</sup>, Ana Bastos<sup>2</sup>, Philippe Ciais<sup>3,4</sup>, Xuhui Wang<sup>1</sup>, Christian Rödenbeck<sup>5</sup>, Pierre Gentine<sup>6</sup>, Frédéric Chevallier<sup>3</sup>, Vincent W. Humphrey<sup>7</sup>, Chris Huntingford<sup>8</sup>, Michael O'Sullivan<sup>9</sup>, Sonia I. Seneviratne<sup>10</sup>, Stephen Sitch<sup>11</sup> & Shilong Piao<sup>1,12✉</sup>

Global fluctuations in annual land carbon uptake (NEE<sub>IAV</sub>) depend on water and temperature variability, yet debate remains about local and seasonal controls of the global dependences. Here, we quantify regional and seasonal contributions to the correlations of globally-averaged NEE<sub>IAV</sub> against terrestrial water storage (TWS) and temperature, and respective uncertainties, using three approaches: atmospheric inversions, process-based vegetation models, and data-driven models. The three approaches agree that the tropics contribute over 63% of the global correlations, but differ on the dominant driver of the global NEE<sub>IAV</sub>, because they disagree on seasonal temperature effects in the Northern Hemisphere (NH, >25°N). In the NH, inversions and process-based models show inter-seasonal compensation of temperature effects, inducing a global TWS dominance supported by observations. Data-driven models show weaker seasonal compensation, thereby estimating a global temperature dominance. We provide a roadmap to fully understand drivers of global NEE<sub>IAV</sub> and discuss their implications for future carbon-climate feedbacks.

<footnote n="1"><sup>1</sup>Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 100871, China. <sup>2</sup>Department of Biogeochemical Integration, Max Planck Institute for Biogeochemistry, Jena 07745, Germany. <sup>3</sup>Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette 91191, France. <sup>4</sup>The Cyprus Institute, 20 Konstantinou Kavafi Street, Nicosia 2121, Cyprus. <sup>5</sup>Department of Biogeochemical Systems, Max Planck Institute for Biogeochemistry, Jena 07745, Germany. <sup>6</sup>Department of Earth and Environmental Engineering, Columbia University, New York, NY 10027, USA. <sup>7</sup>Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA. <sup>8</sup>U.K. Centre for Ecology and Hydrology, Benson Lane, Wallingford OX10 8BB, UK. <sup>9</sup>College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4 4QF, UK. <sup>10</sup>Department of Environmental Systems Science, Institute for Atmospheric and Climate Science, ETH Zurich, Zurich 8092, Switzerland. <sup>11</sup>College of Life and Environmental Sciences, University of Exeter, Exeter EX4 4RJ, UK. <sup>12</sup>State Key Laboratory of Tibetan Plateau Earth System and Resources Environment (TPESRE), Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100085, China. ✉email: slpiao@pku.edu.cn</footnote>
The global land carbon sink offsets about 30% of anthropogenic CO₂ emissions every year, slowing global warming¹,². However, this sink exhibits substantial inter-annual variability (IAV), which causes the CO₂ growth rate (CGR) to fluctuate³,⁴. Understanding the IAV of the net land carbon sink, hereafter the IAV of net ecosystem exchange (NEE<sub>IAV</sub>), is critical for estimating the future evolution of atmospheric CO₂ with the emergence of carbon-climate feedbacks⁵⁻⁷. Such an understanding is especially useful if climate-driven fluctuations of land carbon sink become more frequent or more intense under climate change, which may impact the capability of the land to offset CO₂ emissions. The dominant drivers of the global NEE<sub>IAV</sub> appear to be inconsistent across studies⁸⁻¹⁰. Some studies have shown that the global NEE<sub>IAV</sub> represented by the IAV of the CGR, is strongly related to the variations of the mean annual temperature¹¹,¹². Other studies underlined the dominance of terrestrial water storage (TWS) and thus water availability for vegetation, even when the temperature (T) was accounted for as a co-predictor⁹,¹³.

This apparent controversy among research findings arises from different assessments of ecosystem responses to water availability or T in different regions and seasons, which combine in complex ways to produce the emergence of global drivers of NEE<sub>IAV</sub>¹⁰,¹⁴⁻¹⁶. Hence depending on the datasets used, their aggregation can lead to differing conclusions about the dominant drivers for the global NEE<sub>IAV</sub>. How the distinct contributions of each region sum to the global signal has not been investigated in detail. Most studies have focused on the response of tropical net ecosystem exchange (NEE) to T and water availability⁴,¹²,¹⁷⁻¹⁹. Researchers have placed far less attention on the contribution of northern ecosystems to the global response, despite evidence identifying large NEE<sub>IAV</sub> over North America and Europe in response to interannual variability of both T and water availability²⁰⁻²³. Furthermore, a full understanding of seasonally varying influences of water availability and T on NEE<sub>IAV</sub> for tropical and northern ecosystems is still lacking¹⁵,²¹,²⁴. Water availability has been suggested to be positively correlated with carbon sink anomalies in most seasons over North America, while seasonally opposite (i.e., compensating) influences of T on NEE may weaken the overall impact of T on annual NEE²⁰. Seasonal compensation effects of T on NEE were reported during recent extreme events in northern regions²⁰,²⁵,²⁶. For example, higher net carbon uptake induced by the warmer T in spring was found to be canceled out by less carbon uptake in the following warmer summer, given a legacy drying of soils from spring to summer²⁵,²⁷. It remains unclear to what extent such a seasonal compensation effect of T on NEE influences the overall global response of NEE<sub>IAV</sub> to annual T anomalies.

Here we assessed how the impacts of TWS and T on the NEE<sub>IAV</sub> vary in different regions and seasons, and how these impacts aggregate in space and time to influence the global annual relationship of NEE with TWS and T. To do so, we used three spatially and monthly explicit estimates of NEE. These datasets are (1) three atmospheric inversions (CAMS, Jena CarboScope, and a new Jena CarboScope inversion called NEE-T-WTS)²⁸⁻³¹, (2) dynamic global vegetation models (DGVMs) of the TRENDY project version 7 (refs. ¹,³²), and (3) three machine learning models upscaling NEE from in situ eddy-covariance measurements (FLUXCOM)¹⁰,³³. All of these datasets cover the period from 1979 to 2016. Although the NEE from atmospheric inversions includes emissions from disturbances such as fires, which are not included in FLUXCOM models³³ and some TRENDY models¹, fire emissions are known to have a relatively small impact on the global NEE<sub>IAV</sub>⁴. We focused on determining the contribution from tropical and northern ecosystems to the dependence of global NEE<sub>IAV</sub> on TWS and T. For tropical ecosystems, we examined the contributions of NEE<sub>IAV</sub> during the dry and the wet seasons. For northern ecosystems, we considered four

seasons separately (boreal spring, summer, autumn, and winter). Throughout this paper, we used the sign convention that a positive value of NEE is for a net carbon release from land to the atmosphere.

## Results

Annual global NEE<sub>IAV</sub> correlations with TWS and T (r<sub>TWS</sub> and r<sub>T</sub>). We first calculated the correlation between the annual global NEE<sub>IAV</sub> and the IAV of TWS (r<sub>TWS</sub>) or T (r<sub>T</sub>), from CGR observations (see Methods) and three different NEE datasets listed above, using reconstructed TWS based on the Gravity Recovery and Climate Experiment (GRACE) satellite observations³⁴ and T from a gridded record of meteorological measurements³⁵. We found that at the annual timescale, the global NEE<sub>IAV</sub> was more strongly correlated with TWS than with T when using data from CGR observations, atmospheric inversions, and DGVMs from TRENDY project (Fig. 1a). This result agrees with the findings of Humphrey et al.¹³, which were based on CGR time series as a proxy of the global NEE<sub>IAV</sub>. The strong correlation of NEE with TWS may, however, contain some response to T because of the strong covariation between TWS and T. These covariations may occur especially over shorter timescales because drought and heatwaves usually co-occur as compound events³⁶. Furthermore, atmospheric feedbacks induced by soil moisture deficits produce increases in T and atmospheric dryness that affect the global NEE<sub>IAV</sub>⁹. To filter out the impact of T on the relationship between NEE and TWS and vice versa, we calculated the partial correlations between the global NEE<sub>IAV</sub> and the individual IAVs of T and TWS, after controlling for the effect of the other variable. These results also show that the global NEE<sub>IAV</sub>, as estimated by CGR

<figure><img src="image_4.png" /></figure>

**Fig. 1** Relationships between the interannual variability (IAV) of global net ecosystem exchange (NEE) and terrestrial water storage (TWS) or temperature (T). The correlation (a) and partial correlation (b) between the IAV of global NEE and the IAVs of global mean TWS (r<sub>TWS</sub>) or T (r<sub>T</sub>) are estimated. The IAV of CO₂ growth rate (CGR) is a proxy of the IAV of global NEE. The NEE is also estimated by atmospheric inversions, dynamic global vegetation models (DGVMs) from the TRENDY project in simulation S2 (NEE output; see Methods) and in simulation S3 (NBP output with the opposite sign; see Methods), and FLUXCOM models. Cyan (red) circles indicate the correlations between global NEE and TWS (T). Error bars for atmospheric inversions and FLUXCOM models indicate the range of models, while the error bars for DGVMs from the TRENDY project are the 1-σ inter-model spread. The NEE derived from FLUXCOM models and DGVMs from the TRENDY project in S2 (NEE output) do not include disturbances (e.g., fire emissions), but the NEE in atmospheric inversions and DGVMs from the TRENDY project in S3 (NBP output with the opposite sign) include disturbances. A larger NEE indicates less net carbon uptake.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Contribution of three latitudinal bands to the relationships between global interannual variability of net ecosystem exchange (NEE<sub>IAV</sub>) and terrestrial water storage ($r_{TWS}$) or temperature ($r_T$). In panel a, positive values mean that the NEE<sub>IAV</sub> in the latitudinal band supports the dominance of terrestrial water storage in the NEE<sub>IAV</sub> while the negative values mean that the NEE<sub>IAV</sub> in the latitudinal band supports the dominance of temperature. In panel b, the solid bars are the contributions to the global $r_{TWS}$ ($C^{TWS}$) while open bars are the contributions to the global $r_T$ ($C^T$). Positive values of $C^{TWS}$ or $C^T$ mean that removing the regional NEE<sub>IAV</sub> would reduce the absolute value of the global $r_{TWS}$ or $r_T$. In panels a and b, green bars indicate NEE<sub>IAV</sub> is estimated by atmospheric inversions, while blue bars and brown bars indicate that NEE<sub>IAV</sub> is estimated by DGVMs from the TRENDY project in simulation S2 (NEE output) and by FLUXCOM models, respectively. For atmospheric inversions and FLUXCOM models, the error bars indicate the range of models while for DGVMs from the TRENDY project, the error bars indicate the 1-σ inter-model spread.</figcaption></figure>

observations, atmospheric inversions, and DGVMs, is more strongly linked to TWS than T (Fig. 1b).

By contrast, both simple correlation and partial correlation results show that the global $r_T$ value is larger than the absolute value of global $r_{TWS}$ when using the NEE values of FLUXCOM models (Fig. 1). The more substantial control of temperature on FLUXCOM-estimated global NEE<sub>IAV</sub>, which is opposite to results from atmospheric inversions and DGVMs, links to a higher NEE<sub>IAV</sub> correlation with T rather than a lower correlation with TWS. Indeed, the value of global $r_{TWS}$ obtained from FLUXCOM models ($r_{TWS} = -0.67$) lies in the range of the $r_{TWS}$ values from other datasets in Fig. 1a. In addition, the difference between the simple and partial correlation values with T is generally much larger than with TWS, except for FLUXCOM models (Fig. 1a versus b). This finding suggests a noticeable contribution of TWS to the correlation between global NEE<sub>IAV</sub> and T. This contribution of TWS is also supported by previous results from CGR observations and coupled land-atmosphere models<sup>9,13</sup>.

**Seasonal compensation effects in northern ecosystems.** To further investigate the relative dominance of TWS and T in the overall global NEE<sub>IAV</sub>, we disaggregated the annual global NEE<sub>IAV</sub> correlations with TWS and T ($r_{TWS}$ and $r_T$) into contributions from different regions and seasons (named $C^{TWS}$ and $C^T$). Our disaggregation method exploits the additive properties of covariances (see Methods). Briefly, the sum of all local and seasonal $C^{TWS}$ and $C^T$ values explains the absolute value of the global $r_{TWS}$ and $r_T$, respectively. Positive values of $C^{TWS}$ or $C^T$ mean that removing NEE from a given region/season reduces the absolute value of the global $r_{TWS}$ or $r_T$. As shown in Fig. 2a, all three approaches used to diagnose regional NEE agree on a larger $C^{TWS}$ than $C^T$ in the tropics and the southern extra-tropics. Therefore, the higher global $r_T$ value than the absolute value of $r_{TWS}$ diagnosed from FLUXCOM models must be the consequence of a distinct modulation of $C^{TWS}$ and $C^T$ in the

Northern Hemisphere (NH, >25°N). This leads us to examine more in-depth the contribution of NEE<sub>IAV</sub> in the NH to the global $r_{TWS}$ and $r_T$. We found that FLUXCOM models did estimate a larger $C^T$ for the NH ($C^T = 0.23 \pm 0.01$) than atmospheric inversions ($C^T = 0.12 \pm 0.04$) and DGVMs from TRENDY project in their simulation S2 ($C^T = 0.03 \pm 0.11$) (Fig. 2b). In addition, the $C^T$ of northern ecosystems in FLUXCOM models is also much higher than $C^{TWS}$, while in atmospheric inversions and DGVMs, the $C^T$ of northern ecosystems is almost equal to $C^{TWS}$ (Fig. 2a). This finding clearly illustrates that the dominance of T (rather than TWS) in the global NEE<sub>IAV</sub> estimated by FLUXCOM models takes its origins in the NH.

To get more insights into this NH phenomenon, we then assessed how the large annual $C^T$ of northern ecosystems suggested by FLUXCOM models emerges from different seasons (boreal spring, summer, autumn, and winter). As shown in Fig. 3, all three approaches agree that $C^T$ in summer is the largest contributor to the annual $C^T$ in the NH, as it is much higher than $C^T$ in autumn and winter. This implies a relatively weak role of ecosystem respiration during the non-growing season in the annual NEE<sub>IAV</sub>. The high $C^T$ of the NH in summer is however partially compensated by $C^T$ in spring from all three datasets (Fig. 3). Compared to atmospheric inversions and DGVMs from the TRENDY project in simulation S2, FLUXCOM models instead estimate a much weaker spring-summer compensation of seasonal $C^T$ in northern ecosystems. Specifically, the FLUXCOM-estimated $C^T$ value in summer is larger than in atmospheric inversions and DGVMs, while in spring, it is much weaker with a near-zero value (Fig. 3c). This explains the emerging large annual $C^T$ in northern extra-tropical regions (Fig. 2), and thus indirectly the T-dominance in the global NEE<sub>IAV</sub> in FLUXCOM models by adding to the tropical influence (Fig. 1).

We have identified that the strength of the spring-summer compensation of $C^T$ in northern ecosystems is largely related to the dominant driver of the global NEE<sub>IAV</sub>. Therefore, we next investigated more systematically the magnitude of $C^T$ in spring
<figure><img src="image_6.png" /><figcaption>**Fig. 3** C**T**WS and C**T** in the northern extra-tropics during four seasons. C**T**WS and C**T** indicate contributions to the global correlations of interannual variability of net ecosystem exchange (NEE<sub>IAV</sub>) against terrestrial water storage (r<sub>TWS</sub>) and temperature (r<sub>T</sub>), respectively. NEE<sub>IAV</sub> is estimated by atmospheric inversions (a), DGVMs from the TRENDY project in simulation S2 (NEE output) (b), and FLUXCOM models (c). The four seasons are boreal spring (MAM, March-May), summer (JJA, June-August), autumn (SON, September-November), and winter (DJF, December and January-February). Positive values of C**T**WS or C**T** mean that removing the seasonal NEE<sub>IAV</sub> would reduce the absolute value of the global correlation r<sub>TWS</sub> or r<sub>T</sub>. The larger rectangles indicate the higher absolute values of the C**T**WS or C**T**. For atmospheric inversions and FLUXCOM models, the black edges of the rectangles indicate that the signs of the C**T**WS or C**T** are the same among all models. For DGVMs from the TRENDY project, black edges of the rectangles indicate that the signs of C**T**WS or C**T** derived from more than 10 out of 14 models are consistent with those from the model ensemble mean.</figcaption></figure>

and summer across the NH. To do so, we examined the sensitivities of NEE to T in boreal spring and summer. Results show that in summer, all three approaches agree on a smaller net carbon uptake with warmer T anomalies at northern mid-latitudes (Fig. 4a-c; brown or blue areas), but they are not consistent at northern high latitudes (Fig. 4a-c). On the one hand, atmospheric inversions and DGVMs from TRENDY project in simulation S2 tend to suggest a negative T sensitivity (warmer summer, lower NEE, and more uptake) in the boreal and arctic ecosystems of North Asia (Fig. 4a, b; pink or green areas), which is the opposite of their results over the mid-latitudes. This estimation of negative T sensitivity by DGVMs could be influenced by their lack of permafrost carbon dynamics and nutrient limitation in their simulations<sup>1,37,38</sup>. Atmospheric inversions may be also limited by relatively sparse atmospheric CO<sub>2</sub> measurement stations in northern Eurasia (Supplementary Note 1). Nevertheless, a negative summer T sensitivity of NEE in high-latitude ecosystems is not implausible.

On the other hand, FLUXCOM models indicate that summer T effects on NEE at high latitudes have the same sign as those in the mid-latitude regions (warmer summer, higher NEE, and less uptake; Fig. 4c; brown areas), which could explain the higher overall summer C<sup>T</sup> of northern ecosystems than that estimated by atmospheric inversions and DGVMs from the TRENDY project (Figs. 3c and 4c). To what extent this is because FLUXCOM models are trained by dense flux tower data in temperate regions and extrapolate temperate NEE responses to the high-latitude regions needs to be further studied. A recent study based on eddy-covariance measurements found that higher summer temperature could enhance the net carbon uptake in the Arctic<sup>39</sup>, more consistent with atmospheric inversions and DGVMs rather than FLUXCOM models. Given the diversity of boreal and arctic systems<sup>40-42</sup>, we note that the response of NEE to summer T in the high-latitude regions still remains uncertain.

Compared with summer, the sensitivity of spring NEE to spring T anomalies estimated by FLUXCOM models is different from estimates by the other two approaches almost throughout the NH. In general, atmospheric inversions and DGVMs from the TRENDY project show a negative sensitivity of NEE to spring T (warmer spring, lower NEE, and more uptake) (Fig. 4a, b; green or blue areas), while the T sensitivity from FLUXCOM models is positive (warmer spring, higher NEE, and less uptake) (Fig. 4c; brown areas). Such a widespread positive sensitivity of NEE to spring T estimated by FLUXCOM models is surprising because a warmer spring (i.e., higher T) is commonly observed to enhance

the net carbon sink at most long-term eddy-covariance sites in northern ecosystems<sup>16,43</sup>. This FLUXCOM-estimated positive sensitivity may be due to the temporal extrapolation of the T sensitivity from some short-term sites, where the impact of T on NEE in recent years could be different from the long-term average of the impact of T<sup>44,45</sup>. There may also be an effect whereby the training of FLUXCOM models by both space and time gradients of NEE, may cause subtropical sites to overly affect the prediction of the T sensitivity at temperate and boreal sites<sup>46</sup>. The positive sensitivity of NEE to spring T estimated by FLUXCOM models results in the spring C<sup>T</sup> being positive in some regions of NH, which weakens the overall compensation effects of spring C<sup>T</sup> on summer C<sup>T</sup> (Fig. 3c). Hence FLUXCOM models estimate the global NEE<sub>IAV</sub> to be more correlated with T instead of TWS (Fig. 1).

The results from atmospheric inversions and DGVMs from the TRENDY project demonstrate that the seasonal (spring vs. summer) compensation of the effects of T anomalies on NEE<sub>IAV</sub> in northern ecosystems weakens the control of global T on global NEE<sub>IAV</sub>. This seasonal compensation also weakens the influence of annual T on NEE<sub>IAV</sub> in the NH, thereby explaining the previous findings that the control of global T on global NEE<sub>IAV</sub> is weaker than that of tropical T<sup>4,13</sup>. The seasonal compensation in NH systems occurs mainly between the boreal spring and summer, as indicated by atmospheric inversions and DGVMs from the TRENDY project (Fig. 3). A warmer spring can stimulate carbon uptake in the NH<sup>16,47,48</sup>, causing the contribution of the NH in spring to the global r<sub>T</sub> to be negative. In contrast, warmer summers with potentially higher air temperature than the optimum of vegetation productivity or reduced soil moisture can suppress the net carbon uptake in northern ecosystems<sup>49-52</sup>. The reduced net carbon uptake (higher NEE) in warmer summers results in a positive contribution of NEE<sub>IAV</sub> in northern ecosystems to the global r<sub>T</sub>. In addition, atmospheric inversions and DGVMs from the TRENDY project indicate that the seasonal compensation of the effects of T on NEE<sub>IAV</sub> is not limited to local extreme events such as severe droughts<sup>20,25,27</sup>, but is also a widespread feature arising during normal years in the NH. This compensation of T effects within the NH explains why at the global scale, TWS is the strongest driver of annual variations in NEE, except in FLUXCOM models (Fig. 1). The possible underestimation of NH seasonal compensation effects in FLUXCOM models leads to a dominant role of T rather than TWS in the global NEE<sub>IAV</sub> from this approach, which is in contrast to the results obtained using the observed CGR<sup>13</sup> (Fig. 1).
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Sensitivity of net ecosystem exchange (NEE) to temperature (T) and terrestrial water storage (TWS) in boreal spring and summer. NEE is estimated by atmospheric inversions (a, d), DGVMs from the TRENDY project in simulation S2 (NEE output) (b, e), and FLUXCOM models (c, f). A positive sensitivity to TWS means drier conditions are correlated with more net carbon uptake, while a positive sensitivity to T means warmer conditions are correlated with less net carbon uptake.</figcaption></figure>

**Contribution of tropical NEE during the dry and wet seasons.**

Although the contribution of NEE in northern ecosystems to the global rT plays a critical role in determining the dominant driver of the global NEEIAV, the magnitude of both the global rTWS and rT is mainly contributed by the tropical ecosystems (71–90% and 63–90%, respectively; Fig. 2b). Interannual variability of NEE in tropical ecosystems is much larger than that in the extra-tropics, and thus contributes more to the global NEEIAV²⁻⁴. Given that the impacts of climate variations on carbon dynamics in the tropics depend on hydrothermal conditions⁸, to gain more understanding of annual tropical CTWS and CT, we separated these two contributions into dry and wet seasons (see Methods). Results from all three approaches show that, overall, tropical CTWS has similar value between the dry and the wet seasons, and similarly, tropical CT has similar value between the dry and the wet seasons (Supplementary Fig. 1). Moreover, in each tropical continent, albeit with uncertainties, CTWS and CT during the dry season are also similar to those during the wet season, respectively (Supplementary Fig. 2). Further, the spatial distributions of the absolute values of tropical CTWS and CT during the dry and wet seasons are strongly related to the magnitude of NEEIAV (Fig. 5 and Supplementary Fig. 3). Regions with large absolute values of CTWS and CT shown in Fig. 5 also have large NEEIAV in the same

season. There is a concurrence of large CTWS and strong inter-annual variations of NEE in savannas of southern Africa during the dry season (Fig. 5a–c and Supplementary Fig. 3). This region seems to be a hotspot of NEEIAV, and the large CTWS and NEEIAV during the dry season reflect a tight relationship between NEE and TWS at interannual timescales (Supplementary Fig. 4). More water availability during the dry season can induce a larger enhancement of vegetation productivity over ecosystem respiration, thus facilitating more net carbon uptake in this region²⁴.

We investigated whether the general length of the dry season in tropical regions could partly explain the spatial and seasonal variations of NEEIAV (Supplementary Figs. 3 and 5) and thus those of CTWS and CT. Specifically, we found that regions with large absolute values of CTWS and CT in the dry season have longer dry seasons and similarly for the wet season (Fig. 5 and Supplementary Fig. 5). This finding indicates that future changes in the length of the dry season in tropical regions could affect the spatiotemporal patterns of CTWS and CT (ref. ⁵³). In addition, we found negative values of CTWS or CT in some locations, e.g., for tropical Asia (Fig. 5), which result from the covariations between NEEIAV in those regions and the IAVs of global TWS or global T (see Methods). Two mutually exclusive mechanisms could explain the negative values of CTWS or CT (see Methods). One
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Spatial distribution of the $C^{\text{TWS}}$ and $C^T$ in the tropics during the dry and wet seasons. $C^{\text{TWS}}$ and $C^T$ indicate contributions to the global correlations of interannual variability of net ecosystem exchange (NEE$_{IAV}$) against terrestrial water storage and temperature, respectively. NEE$_{IAV}$ is estimated by atmospheric inversions (a, d), DGVMs from the TRENDY project in simulation S2 (NEE output) (b, e), and FLUXCOM models (c, f). Note that pixels where there is only a dry or only a wet season within a year are excluded. Triangles indicate that using simulated soil moisture by DGVMs or water availability index (WAI) for FLUXCOM models as terrestrial water storage instead of using observation-based terrestrial water storage would change the sign of $C^{\text{TWS}}$.</figcaption></figure>

is that years with more water availability or lower temperature do not enhance carbon sink in those regions. A second possibility is that local water availability or temperature variations are in the opposite direction of the global fluctuations (see Methods). The negative local $C^{\text{TWS}}$ and $C^T$ values imply a spatial compensation with positive values elsewhere in the tropics (Supplementary Note 2), modulating the $C^{\text{TWS}}$ and $C^T$ in the whole tropics.

Discussion

We have shown that the relative magnitude of $C^{\text{TWS}}$ and $C^T$ in the NH is critical to determine the dominant driver of the global NEE$_{IAV}$ (Fig. 2a). However, there are considerable uncertainties in the estimate of $C^{\text{TWS}}$ and $C^T$ in subregions of NH (Supplementary Note 1). Specifically, we found little consensus between the three approaches in estimating $C^{\text{TWS}}$ and $C^T$ for many subregions in the NH (Fig. 6). Such regional discrepancies in annual $C^{\text{TWS}}$ and $C^T$ could partly result from seasonal compensation differences. For example, atmospheric inversions estimate larger $C^{\text{TWS}}$ values than the other two approaches in Asian regions (i.e., North Asia, northern East Asia, Central Asia, South Asia, and subtropical China; Fig. 6), inducing a larger $C^{\text{TWS}}$ in the NH (Fig. 2b). The large annual $C^{\text{TWS}}$ in Asian regions, as estimated

by atmospheric inversions, links to the positive values of spring $C^{\text{TWS}}$ in these regions, which add to $C^{\text{TWS}}$ in other seasons (Supplementary Fig. 6). By contrast, DGVMs from TRENDY project in simulation S2 and FLUXCOM models suggest the negative values of spring $C^{\text{TWS}}$ in Asian regions, weakening the annual $C^{\text{TWS}}$ (Supplementary Fig. 6). DGVMs from TRENDY project and FLUXCOM models suggest that increased water availability in spring can enhance the net land carbon uptake in East Asia and the arctic region of Asia (Fig. 4e, f; green and blue areas), by promoting the photosynthesis of vegetation$^{10,54}$. Yet this is not supported by NEE values derived from atmospheric inversions (Fig. 4d; brown and pink areas).

In the tropics, discrepancies among the three approaches can be seen when analyzing continental contributions to the global $r_{\text{TWS}}$ and $r_T$ (Fig. 6 and Supplementary Fig. 2). Results from the three approaches show that both tropical Africa and tropical America generally contribute more to the global $r_{\text{TWS}}$ and $r_T$ than tropical Asia (including northern Australia, as shown in Fig. 6a). But large uncertainties remain in quantifying the continental contributions of tropical Africa and tropical America (Fig. 6 and Supplementary Fig. 2). Palmer et al.$^{17}$, using inversions based on remotely-sensed column CO$_2$ mole fraction retrievals, found some evidence of a substantial role of tropical Africa in the
<figure><img src="image_9.png" /><figcaption>**Fig. 6** C**TWS** and C**T** in tropical continents and northern extra-tropical subregions. C**TWS** and C**T** indicate contributions to the global correlations of interannual variability of net ecosystem exchange against terrestrial water storage and temperature, respectively. As shown in panel **a**, the tropics is divided into three continents (Trop.Am tropical America, Trop.Af tropical Africa, Trop.As tropical Asia, with the latter including northern Australia). Northern Hemisphere (NH, >25°N) is divided into ten subregions. BAm Boreal America, NEu North Europe, NAs North Asia, TEu Temperate Europe, WAm Western North America, EAm Eastern North America, Med Mediterranean, CAs Central Asia, SAs South Asia and subtropical China, and NEAs Northern East Asia. For atmospheric inversions and FLUXCOM models, the black edges of the rectangles in panel **b** indicate that the signs of the C**TWS** or C**T** are the same among all models. For DGVMs from the TRENDY project, the black edges of the rectangles in panel **b** indicate that the signs of C**TWS** or C**T** derived from more than 10 out of 14 models are consistent with those from the model ensemble mean.</figcaption></figure>

variability of the overall tropical NEE. However, that study was limited to a period of four years only. The longer NEE time series of atmospheric inversions based on in situ surface network used here show the large contributions of tropical Africa hold for the past four decades. In tropical Africa, fire is not likely to be the cause of net carbon flux anomalies in the recent few years¹⁷,⁵⁵, implying that the control of NEEIAV is the balance between photosynthesis and respiration. Besides the large contributions of tropical Africa, DGVMs and FLUXCOM models also estimate large C**TWS** and C**T** of tropical America (Fig. 6). This is probably because of the high NEEIAV in tropical America estimated by DGVMs and FLUXCOM models⁴,¹⁰, reflecting the high sensitivity/vulnerability of tropical forests to climate anomalies⁵⁶. The uncertainty in the contributions of three tropical continents may be largely from the spatial compensation effects (Supplementary Note 2).

To provide more accurate assessments of the spatial compensation for C**TWS**, or for C**T**, a key priority is reducing the uncertainties of the three approaches in estimating tropical NEEIAV. Atmospheric inversions assimilating surface measurements of atmospheric CO₂ have potentially low capacity to separate the IAV signal between tropical continents, given the sparsity of stations recording atmospheric CO₂ in the tropics⁵⁷. The NEE estimate by the three atmospheric inversions is thus likely influenced by the prior settings in tropical regions. This can be partly remedied by assimilating total CO₂ dry air mole fraction

retrievals from satellites, as these have better coverage of the tropics³⁰,⁵⁸. For the FLUXCOM models, besides limitations from relatively few training samples in the tropics, the estimation of NEE is affected by both specific machine learning methods and in particular, the calculation of the water availability index (WAI) used as a predictor of NEE³³,⁴⁶. The soil water storage capacity in the bucket model used to simulate WAI in FLUXCOM is constant and relatively small⁴⁶. By contrast, DGVMs calculate their own soil moisture, aiming to better capture the actual water availability and water stress on NEE. Uncertainties in simulating the NEE variability in DGVMs are thus mainly derived from the model structure, including soil hydrology and limitations of soil moisture and atmospheric dryness on plant photosynthesis and soil respiration³²,⁵⁹,⁶⁰.

Atmospheric inversions and DGVMs from the TRENDY project indicate a dominant role of TWS on the global NEEIAV. However, TWS anomalies in some tropical regions relate to water storage in lakes, floodplains, and wetlands, which should have less effect on the carbon uptake than soil moisture⁶¹. This suggests a stronger relationship of NEEIAV with simulated soil moisture in DGVMs or with WAI in FLUXCOM models than with observation-based TWS in tropical regions (Supplementary Figs. 7 and 8). This finding is also because these two model-based approaches are more consistent with their own derived soil moisture than with the independently observed TWS. Nevertheless, results from DGVMs and FLUXCOM models when considering their own
simulated soil moisture, also support that the dry season length controls the absolute values of C**TWS** in the tropics (Fig. 5 and Supplementary Fig. 9), although showing different spatial compensation patterns compared to results where C**TWS** is calculated from observation-based TWS (Fig. 5; Supplementary Fig. 9; and Supplementary Note 2). In addition, even though the correlation between IAV of TWS reconstructed from GRACE satellite observations and IAV of the simulated soil moisture or water availability by the models is relatively weak in boreal and Arctic regions of the NH (Supplementary Fig. 10), it has little influence on the assessment of the dominant driver of NEE IAV in the NH (Supplementary Fig. 7).

Our study revealed how the interannual NEE variations in different regions and seasons contribute to the overall apparent controls of TWS and T on the global NEE. At the global scale, both atmospheric inversions and DGVMs suggest that TWS exerts a stronger control on the NEE than T. This finding of a stronger TWS-dependence is supported by the tight relationship between interannual variations of the data-led global NEE estimated by atmospheric inversions and observed CGR variations⁴. The weaker T control primarily results from the seasonal compensation that the negative impacts of T on spring NEE offset the positive impacts on summer NEE over the northern ecosystems. However, a weaker spring-summer compensation of T effects in FLUXCOM models instead makes them produce a stronger T control than TWS on the global NEE. We note that the weak seasonal compensation is not supported by local flux tower measurements, which show a larger net carbon uptake in land in a warmer spring across northern ecosystems¹⁶,⁴³. Nevertheless, weaker seasonal compensation effects over northern ecosystems may emerge in the future, due to the reduced chilling or growing light limitations in spring⁴⁴ and to the accelerating carbon release from permafrost thawing or from more frequent summer heatwaves³⁸,⁶². A weaker seasonal compensation effect in northern ecosystems could strengthen the T control in the global land carbon uptake variability. Hence monitoring how the seasonal compensation effects are changing in the Northern Hemisphere, supported by improved process modeling, will help better inform society about the emergence of any altered future carbon-climate feedbacks.

Atmospheric inversions, DGVMs from the TRENDY project, and FLUXCOM models all agree on the importance of the tropical contributions to the global NEE-TWS and NEE-T relationships. Here we show that this pan-tropical influence results from almost equal contributions from the dry and wet seasons. Our study highlights the role of dry season length in determining the spatial and seasonal pattern of the tropical contribution. To eventually quantify the contribution of tropical NEE during different seasons and at a finer spatial scale, a more accurate assessment of the spatial compensation of the contributions is needed. At present, the accuracy of tropical NEE estimates from atmospheric inversions and FLUXCOM models at fine resolution is limited by sparse atmospheric CO₂ observations and eddy-covariance measurements in the tropics, respectively³³,⁵⁷. Assimilating satellite measurements of total CO₂ dry air mole fraction into inversion algorithms has the potential to improve NEE estimates of atmospheric inversions at small spatial scales³⁰,⁵⁸. For DGVMs, the poor calibration or even lack of representation of key carbon cycle processes, such as nutrient limitations and access of roots to groundwater, could potentially induce biases in estimating NEE variations³⁷,⁶³. Given these different sources of uncertainty in simulating NEE among the three approaches, we emphasize the importance of comparing independent data sources and methods in determining regional carbon cycle attributes, as done in this study. In summary, we have shown that global annual variations in NEE may be more linked

to fluctuations in terrestrial water storage than temperature, although this depends on a selected strand of evidence. Hence, we determined where data streams and model projections differ in estimating seasonal and regional NEE, and offsets in their variations. Our research provides guidance for future measurement campaigns and model development, which will lower further uncertainty on the dominant drivers of large-scale variations in NEE.

## Methods

**Atmospheric CO₂ inversions.** We used gridded net carbon fluxes for the period of 1979–2016 from two long-term atmospheric inversion models: CAMS version 17r1 (ref. ²⁸) and Jena CarboScope version s76oc_v2020 (ref. ²⁹). The two inversions both assimilated surface measurements of atmospheric CO₂. The spatial resolutions of net carbon flux from CAMS and Jena CarboScope are 1.9°latitude × 3.75°longitude and 4°latitude × 5°longitude, respectively. We used the remapped monthly output with a spatial resolution of 1° × 1°. We also used another Jena CarboScope inversion called NEE-T-TWS (version sEXTocNEETWS_79r18_v2020) for the period of 1979–2016. In the NEE-T-TWS inversion, a multilinear regression term related to T and TWS with inversely adjusted coefficients replaces the interannual NEE variations¹⁶,³¹. For some comparison, we used the CAMS inversion assimilating total column CO₂ dry air mole fraction retrievals from the Japanese Greenhouse gases Observing SATellite (GOSAT) for the time period of 2010–2016 (ref. ³⁰). For comparison with the inversion based on GOSAT retrievals, CAMS assimilating surface measurements is also limited to 2010–2016 (named Surf). We defined the net carbon flux as NEE and a positive value of NEE as a net carbon release from land.

**Dynamic global vegetation models.** We used monthly net carbon fluxes estimated by 14 dynamic global vegetation models (DGVMs) from TRENDYv7 project¹,³². The 14 models are CABLE-POP⁶⁴, CLASS-CTEM⁶⁵, CLM5.0⁶⁶, DLEM⁶⁷, ISAM⁶⁸, JSBACH⁶⁹, LPJ⁷⁰, LPX-Bern⁷¹, OCN⁷², ORCHIDEE-CNP⁷³, ORCHIDEE-Trunk⁷⁴, SDGVM⁷⁵, SURFEX⁷⁶, and VISIT⁷⁷, which have a monthly output covering the period of 1979–2016 (ref. ¹). Nine of the 14 models include a fire module. The effects of carbon-nitrogen interactions are considered in nine models while few models include the phosphorus-cycle module such as ORCHIDEE-CNP. There is no model in TRENDYv7 specifically considering the tropical wetland hydrology. For permafrost, no TRENDYv7 DGVMs have explicitly considered it in the simulations. In addition, each model includes forestation activities and afforestation activities or forest regrowth after cropland abandonment, as well as more generally land cover change. For a complete summary of the models, we refer to Le Quéré et al.¹. Three simulations (S1, S2, and S3) were designed in the TRENDY protocol. DGVMs in S3 are forced by varying CO₂, climate change, and land cover change while those in S2 are forced by varying CO₂ and climate change. Differences of carbon fluxes between S2 and S3 are thus from the effects of land use change. We used the NEE outputs in simulation S2 of the TRENDY protocol and compared the results with the net biome productivity (NBP) output in simulation S3. All variables are remapped to a spatial resolution of 0.5° × 0.5°. We mainly used the IAV of modeled NEE in a comparison with that of the NEE from atmospheric inversions and FLUXCOM models.

**FLUXCOM global carbon flux dataset.** We used data-driven NEE retrievals from the FLUXCOM model ensemble from an upscaling flux tower dataset using satellite and meteorological/climate forcing (RS+METEO)¹⁰,³³,⁴⁶. The meteorological forcing dataset used is CRU-NCEP v8 (ref. ⁷⁸). FLUXCOM models use several products from the Moderate Resolution Imaging Spectroradiometer (MODIS), including MOD11A2 land surface temperature, the MOD13A2 vegetation index, MOD15A2 leaf area index and fraction of absorbed photosynthetic active radiation, and bidirectional reflectance distribution function-corrected surface reflectance of MCD43A2 and MCD43A4 (ref. ⁴⁶), which can be obtained from http://daac.ornl.gov/MODIS/. The water availability index (WAI), which is calculated by a bucket model³³,⁴⁶, is also used in driving the machine learning algorithms. The bucket model considers a simple water budget without changes in snow cover, permafrost, and runoff¹,³,⁴⁶. FLUXCOM estimates gridded NEE based on satellite and meteorological/climate forcing data, using three machine learning methods (random forests, artificial neural networks, and multivariate adaptive regression splines) to upscale the NEE from measurements at eddy-covariance sites⁴⁶. The machine learning is based on both the spatial and temporal gradients of observed NEE. More accurate predictions of NEE from FLUXCOM models can be obtained in regions with more eddy-covariance sites, such as temperate and boreal forests⁴⁶. More observations in the tropical and arctic ecosystems are needed. We used the monthly NEE output of FLUXCOM models with a spatial resolution of 1° × 1° and covering the period of 1979–2016.

**Climate and terrestrial water storage data.** We mainly used terrestrial water storage (TWS) as a proxy for water availability. The TWS dataset used in this study is reconstructed using a statistical model based on GRACE observations of TWS
(GSFC dataset), and the relationship between climate drivers (MSWEP precipita- tion dataset and ERA5 air temperature dataset) and TWS³⁴. Monthly reconstructed TWS covers the period of 1979–2016 with a spatial resolution of 1° × 1°. The main analyses in this study (Figs. 1–6) are based on this reconstructed TWS time series. For some comparisons, we also used the simulated soil moisture by DGVMs and the water availability index in FLUXCOM models derived from a soil water balance model forced by precipitation and evapotranspiration¹⁰,⁴⁶. The correlation between TWS and simulated soil moisture or water availability index was estimated (Sup- plementary Fig. 10). Monthly gridded air temperature, precipitation, and potential evapotranspiration (PET) data for 1979–2016 are obtained from the Climatic Research Unit (CRU TS4.02), with a spatial resolution of 0.5° × 0.5° (ref. ³⁵).

**Dry and wet seasons.** The dry season in each grid cell is defined as the period when monthly PET is larger than precipitation in this pixel. The period in which monthly PET is lower than precipitation is defined as the wet season. Our analysis of the impacts of TWS and T on the tropical NEEIAV during the dry and wet seasons excluded the pixels where there is only dry or only wet season within a year (Supplementary Fig. 5).

**CO₂ growth rate (CGR).** We used atmospheric CO₂ mole fraction data from the Greenhouse Gas Marine Boundary Layer Reference of the National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA/ESRL)⁷⁹. This dataset is constructed using measurements of weekly air samples from the Cooperative Air Sampling Network. We calculated the difference of CO₂ mole fraction in January of this year and the following year to obtain the annual CGR.

**Correlation between global annual NEE and TWS or T (r<sub>TWS</sub> or r<sub>T</sub>).** To obtain the IAV of TWS, T, and NEE, we removed the long-term trend using linear regression from the monthly time series and summed the values of each year into the yearly anomalies. Then the Pearson's correlation coefficient between yearly NEE and yearly TWS (r<sub>TWS</sub>) or T (r<sub>T</sub>) was calculated. We excluded the years 1982, 1991, and 1992 from the time series of NEE, T, and TWS for each pixel because volcanic eruptions affected the radiation and the relationship between carbon uptake and water availability.

**Contributions of regional and seasonal NEE to the global r<sub>TWS</sub> and r<sub>T</sub>.** We first analyzed the contribution of monthly NEE in each grid cell to the global r<sub>TWS</sub> and r<sub>T</sub>. We considered extensive quantities for NEE so we can split the global annual NEE anomaly (X<sub>G</sub>, unit: PgC yr⁻¹) into n grid cells and 12 months covering the entire land surface and whole years with Σ<sup>n</sup><sub>i=1</sub>Σ<sup>12</sup><sub>m=1</sub> x<sub>i,m</sub> = X<sub>G</sub>, where the x<sub>i,m</sub> is the NEE anomaly (unit: PgC month⁻¹) in grid cell i and month m. The correlation coefficient between X<sub>G</sub> and global mean T anomaly (T<sub>G</sub>) can be decomposed as:

$$ \begin{aligned} \operatorname{corr}(X_G, T_G) &= \frac{\operatorname{cov}(X_G, T_G)}{\sigma_{X_G, T_G}} = \sum_{i=1}^{n} \frac{\sum_{m=1}^{12} \operatorname{cov}(x_{i,m}, T_G)}{\sigma_{X_G, T_G}} \\ &= \sum_{i=1}^{n} \sum_{m=1}^{12} \frac{\operatorname{corr}(x_{i,m}, T_G) \sigma_{x_{i,m}}}{\sigma_{X_G, T_G}} = \sum_{i=1}^{n} \sum_{m=1}^{12} \operatorname{corr}(x_{i,m}, T_G) \frac{\sigma_{x_{i,m}}}{\sigma_{X_G}} \end{aligned} \quad (1) $$

where cov(X<sub>G</sub>, T<sub>G</sub>) is the covariance between X<sub>G</sub> and T<sub>G</sub>. The σ<sub>X<sub>G</sub></sub> and σ<sub>T<sub>G</sub></sub> are the standard deviation of X<sub>G</sub> and T<sub>G</sub>, respectively. The global r<sub>T</sub> is thus the sum of the correlations between x<sub>i,m</sub> and T<sub>G</sub> weighted by the ratio between the standard deviation of x<sub>i,m</sub> and the standard deviation of X<sub>G</sub> (Eq. (1)). The contribution of x<sub>i,m</sub> to the global r<sub>T</sub> can be expressed as:

$$ C_{i,m}^{T} = \operatorname{corr}(x_{i,m}, T_G) \cdot \frac{\sigma_{x_{i,m}}}{\sigma_{X_G}} \quad (2) $$

which relates to the global correlation as:

$$ \operatorname{corr}(X_G, T_G) = \sum_{i=1}^{n} \sum_{m=1}^{12} C_{i,m}^{T} \quad (3) $$

The method of calculating the contribution of x<sub>i,m</sub> to the global r<sub>TWS</sub> is similar but with a negative sign because the global r<sub>TWS</sub> is negative.

$$ C_{i,m}^{TWS} = -\operatorname{corr}(x_{i,m}, TWS_G) \cdot \frac{\sigma_{x_{i,m}}}{\sigma_{X_G}} \quad (4) $$

where TWS<sub>G</sub> is the global mean TWS anomaly. The magnitudes of C<sup>TWS</sup><sub>i,m</sub> and C<sup>T</sup><sub>i,m</sub> are mainly affected by the ratio between the standard deviation of x<sub>i,m</sub> and the standard deviation of X<sub>G</sub>, especially in the tropics. The magnitudes of the correlation between local NEEIAV and local TWS or T have a weak influence on the spatiotemporal patterns of C<sup>TWS</sup> or C<sup>T</sup> in tropical ecosystems (Fig. 5, Supplementary Figs. 4 and 11). The sign of C<sup>TWS</sup><sub>i,m</sub> or C<sup>T</sup><sub>i,m</sub> is influenced by the correlation between the IAV of local NEE and global TWS or global T, which is probably linked to the control of local TWS or local T in the local NEEIAV, as well as the relationship between the IAVs of climate drivers at local and global scales.

The contribution of annual or seasonal NEE to the global r<sub>TWS</sub> or r<sub>T</sub> can be calculated by the direct sum of monthly C<sup>TWS</sup> or C<sup>T</sup>. For tropical ecosystems, we examined C<sup>TWS</sup> and C<sup>T</sup> during the dry and wet seasons. For northern ecosystems, we considered four seasons: boreal spring (March–May), summer (June–August),

autumn (September–November), and winter (December, January, and February). The contribution of regional or latitudinal NEE to the global r<sub>TWS</sub> or r<sub>T</sub> is calculated by the direct sum of C<sup>TWS</sup> or C<sup>T</sup> in each grid cell of the region or the latitudinal band.

**Sensitivity of NEEIAV to TWS and T for different seasons.** To investigate the sensitivity of NEEIAV to TWS and T in different seasons, we first separated the monthly time series of NEE, TWS, and T into 12 groups, corresponding to the 12 months of a year. Each group included the values for the same month for all years. Then the linear trend of the time series in each group was removed to obtain anomalies of NEE, TWS, and T. Groups of months belonging to the same season were summed to become one group representing the season. The new group thus includes the values of the same season for all years. The sensitivity of NEEIAV to TWS (a<sup>TWS</sup>) and T (a<sup>T</sup>) for a given season (s) and grid cell (i) was estimated by linear multiple regression:

$$ NEE_{i,s} = a_{i,s}^{TWS} \cdot TWS_{i,s} + a_{i,s}^{T} \cdot T_{i,s} + \varepsilon \quad (5) $$

where NEE<sub>i,s</sub>, TWS<sub>i,s</sub>, and T<sub>i,s</sub> are vectors of NEE, TWS, and T in a given season (s) and grid cell (i) during 1979–2016, respectively and ε is a residual term. a<sup>TWS</sup><sub>i,s</sub> and a<sup>T</sup><sub>i,s</sub> indicate the sensitivities of NEEIAV to TWS and T for a given season (s) and grid cell (i), respectively.

## Data availability

The net land carbon fluxes of CAMS atmospheric inversions are available at https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-greenhouse-gas-inversion. The Jena CarboScope inversions results are available from the Jena CarboScope website http://www.BGC-Jena.mpg.de/CarboScope/. The simulations from TRENDY DGVMs are available at https://sites.exeter.ac.uk/trendy. The FLUXCOM ensemble of carbon fluxes is available at www.bgc-jena.mpg.de/geodb/projects/Data.php. Monthly gridded air temperature, precipitation, and potential evapotranspiration data from the Climatic Research Unit can be accessed at https://rudata.uea.ac.uk/cru/data/hrg/. The reconstructed TWS dataset based on GRACE observations is accessible at https://doi.org/10.6084/m9.fishshare.7670849. The CO₂ mole fraction data from the Greenhouse Gas Marine Boundary Layer Reference of the NOAA ESRL is available at https://gml.noaa.gov/ccgg/mbl/data.php.

## Code availability

All computer codes for the analysis of the data are available from the corresponding author on reasonable request.

Received: 2 November 2020; Accepted: 30 May 2022; Published online: 16 June 2022

## References

1. Le Quéré, C. et al. Global carbon budget 2018. Earth Syst. Sci. Data 10, 2141–2194 (2018).

2. Friedlingstein, P. et al. Global carbon budget 2020. Earth Syst. Sci. Data 12, 3269–3340 (2020).

3. Rayner, P. J. et al. Interannual variability of the global carbon cycle (1992–2005) inferred by inversion of atmospheric CO₂ and δ¹³CO₂ measurements. Glob. Biogeochem. Cycles 22, 1–12 (2008).

4. Piao, S. et al. Interannual variation of terrestrial carbon cycle: Issues and perspectives. Glob. Chang. Biol. 26, 300–318 (2020).

5. Betts, R. A. et al. A successful prediction of the record CO₂ rise associated with the 2015/2016 El Niño. Philos. Trans. R. Soc. B Biol. Sci. 373, 20170301 (2018).

6. Keeling, C. D., Whorf, T. P., Wahlen, M. & van der Plicht, J. Interannual extremes in the rate of rise of atmospheric carbon dioxide since 1980. Nature 375, 666–670 (1995).

7. Cox, P. M. et al. Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability. Nature 494, 341–344 (2013).

8. Fang, Y. et al. Global land carbon sink response to temperature and precipitation varies with ENSO phase. Environ. Res. Lett. 12, 064007 (2017).

9. Humphrey, V. et al. Soil moisture–atmosphere feedback dominates land carbon uptake variability. Nature 592, 65–69 (2021).

10. Jung, M. et al. Compensatory water effects link yearly global land CO₂ sink changes to temperature. Nature 541, 516–520 (2017).

11. Wang, W. et al. Variations in atmospheric CO₂ growth rates coupled with tropical temperature. Proc. Natl Acad. Sci. USA 110, 13061–13066 (2013).

12. Wang, X. et al. A two-fold increase of carbon cycle sensitivity to tropical temperature variations. Nature 506, 212–215 (2014).

13. Humphrey, V. et al. Sensitivity of atmospheric CO₂ growth rate to observed changes in terrestrial water storage. Nature 560, 628–631 (2018).
14. Marcolla, B., Rödenbeck, C. & Cescatti, A. Patterns and controls of inter-annual variability in the terrestrial carbon budget. *Biogeosciences* **14**, 3815–3829 (2017).

15. Yin, Y. et al. Changes in the response of the northern hemisphere carbon uptake to temperature over the last three decades. *Geophys. Res. Lett.* **45**, 4371–4380 (2018).

16. Rödenbeck, C., Zaehle, S., Keeling, R. & Heimann, M. How does the terrestrial carbon exchange respond to inter-annual climatic variations? A quantification based on atmospheric CO₂ data. *Biogeosciences* **15**, 2481–2498 (2018).

17. Palmer, P. I. et al. Net carbon emissions from African biosphere dominate pan-tropical atmospheric CO₂ signal. *Nat. Commun.* **10**, 1–9 (2018).

18. Fan, L. et al. Satellite-observed pantropical carbon dynamics. *Nat. Plants* **5**, 944–951 (2019).

19. Poulter, B. et al. Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle. *Nature* **509**, 600–603 (2014).

20. Hu, L. et al. Enhanced North American carbon uptake associated with El Niño. *Sci. Adv.* **5**, 1–11 (2019).

21. Liu, Z. et al. Increased high-latitude photosynthetic carbon gain offset by respiration carbon loss during an anomalous warm winter to spring transition. *Glob. Chang. Biol.* **26**, 682–696 (2020).

22. Reichstein, M. et al. Determinants of terrestrial ecosystem carbon balance inferred from European eddy covariance flux sites. *Geophys. Res. Lett.* **34**, 1–5 (2007).

23. Shiga, Y. P. et al. Forests dominate the interannual variability of the North American carbon sink. *Environ. Res. Lett.* **13**, 084015 (2018).

24. Wang, X., Ciais, P., Wang, Y. & Zhu, D. Divergent response of seasonally dry tropical vegetation to climatic variations in dry and wet seasons. *Glob. Chang. Biol.* **24**, 4709–4717 (2018).

25. Wolf, S. et al. Warm spring reduced carbon cycle impact of the 2012 US summer drought. *Proc. Natl Acad. Sci. USA* **113**, 5880–5885 (2016).

26. Liu, J. et al. Detecting drought impact on terrestrial biosphere carbon fluxes over contiguous US with satellite observations. *Environ. Res. Lett.* **13**, 095003 (2018).

27. Bastos, A. et al. Direct and seasonal legacy effects of the 2018 heat wave and drought on European ecosystem productivity. *Sci. Adv.* **6**, eaba2724 (2020).

28. Chevallier, F. et al. Inferring CO₂ sources and sinks from satellite observations: Method and application to TOVS data. *J. Geophys. Res.* **110**, D24309 (2005).

29. Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO₂ flux history 1982-2001 inferred from atmospheric data using a global inversion of atmospheric transport. *Atmos. Chem. Phys.* **3**, 1919–1964 (2003).

30. Chevallier, F. et al. Objective evaluation of surface- and satellite-driven carbon dioxide atmospheric inversions. *Atmos. Chem. Phys.* **19**, 14233–14251 (2019).

31. Rödenbeck, C., Zaehle, S., Keeling, R. & Heimann, M. The European carbon cycle response to heat and drought as seen from atmospheric CO₂ data for 1999–2018. *Philos. Trans. R. Soc. B Biol. Sci.* **375**, 20190506 (2020).

32. Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon dioxide. *Biogeosciences* **12**, 653–679 (2015).

33. Jung, M. et al. Scaling carbon fluxes from eddy covariance sites to globe: Synthesis and evaluation of the FLUXCOM approach. *Biogeosciences* **17**, 1343–1365 (2020).

34. Humphrey, V. & Gudmundsson, L. GRACE-REC: a reconstruction of climate-driven water storage changes over the last century. *Earth Syst. Sci. Data* **11**, 1153–1170 (2019).

35. Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. *Sci. Data* **7**, 109 (2020).

36. Zhou, S., Zhang, Y., Williams, A. P. & Gentine, P. Projected increases in intensity, frequency, and terrestrial carbon costs of compound drought and aridity events. *Sci. Adv.* **5**, 1–9 (2019).

37. Du, E. et al. Global patterns of terrestrial nitrogen and phosphorus limitation. *Nat. Geosci.* **13**, 221–226 (2020).

38. Turetsky, M. R. et al. Permafrost collapse is accelerating carbon release. *Nature* **569**, 32–34 (2019).

39. Li, Z. L. et al. Changes in net ecosystem exchange of CO₂ in Arctic and their relationships with climate change during 2002–2017. *Adv. Clim. Chang. Res.* **12**, 475–481 (2021).

40. Walker, D. A. et al. The Circumpolar Arctic vegetation map. *J. Veg. Sci.* **16**, 267–282 (2005).

41. Virkkala, A. M. et al. Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties. *Glob. Chang. Biol.* **27**, 4040–4059 (2021).

42. Randazzo, N. A. et al. Higher autumn temperatures lead to contrasting CO₂ flux responses in boreal forests versus tundra and shrubland. *Geophys. Res. Lett.* **48**, e2021GL093843 (2021).

43. Richardson, A. D. et al. Influence of spring and autumn phenological transitions on forest ecosystem productivity. *Philos. Trans. R. Soc. B Biol. Sci.* **365**, 3227–3246 (2010).

44. Piao, S. et al. Weakening temperature control on the interannual variations of spring carbon uptake across northern lands. *Nat. Clim. Chang.* **7**, 359–363 (2017).

45. Peñuelas, J. et al. Shifting from a fertilization-dominated to a warming-dominated period. *Nat. Ecol. Evol.* **1**, 1438–1445 (2017).

46. Tramontana, G. et al. Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms. *Biogeosciences* **13**, 4291–4313 (2016).

47. Randerson, J. T., Field, C. B., Fung, I. Y. & Tans, P. P. Increases in early season ecosystem uptake explain recent changes in the seasonal cycle of atmospheric CO₂ at high northern latitudes. *Geophys. Res. Lett.* **26**, 2765–2768 (1999).

48. Black, T. A. et al. Increased carbon sequestration by a boreal deciduous forest in years with a warm spring. *Geophys. Res. Lett.* **27**, 1271–1274 (2000).

49. Wang, T. et al. Emerging negative impact of warming on summer carbon uptake in northern ecosystems. *Nat. Commun.* **9**, 1–7 (2018).

50. Huang, M. et al. Air temperature optima of vegetation productivity across global biomes. *Nat. Ecol. Evol.* **3**, 772–779 (2019).

51. Seneviratne, S. I. et al. Investigating soil moisture-climate interactions in a changing climate: a review. *Earth-Sci. Rev.* **99**, 125–161 (2010).

52. Buermann, W. et al. Recent shift in Eurasian boreal forest greening response may be associated with warmer and drier summers. *Geophys. Res. Lett.* **41**, 1995–2002 (2014).

53. Fu, R. et al. Increased dry-season length over southern Amazonia in recent decades and its implication for future climate projection. *Proc. Natl Acad. Sci. USA* **110**, 18110–18115 (2013).

54. Chen, C. et al. Identifying critical climate periods for vegetation growth in the northern hemisphere. *J. Geophys. Res. Biogeosci.* **123**, 2541–2552 (2018).

55. Gloor, E. et al. Tropical land carbon cycle responses to 2015/16 El Niño as recorded by atmospheric greenhouse gas and remote sensing data. *Philos. Trans. R. Soc. B Biol. Sci.* **373**, 20170302 (2018).

56. Saatchi, S. et al. Detecting vulnerability of humid tropical forests to multiple stressors. *One Earth* **4**, 988–1003 (2021).

57. Peylin, P. et al. Global atmospheric carbon budget: results from an ensemble of atmospheric CO₂ inversions. *Biogeosciences* **10**, 6699–6720 (2013).

58. Liu, J. et al. Carbon monitoring system flux net biosphere exchange 2020 (CMS-Flux NBE 2020). *Earth Syst. Sci. Data* **13**, 299–330 (2021).

59. Quetin, G. R., Bloom, A. A., Bowman, K. W. & Konings, A. G. Carbon flux variability from a relatively simple ecosystem model with assimilated data is consistent with terrestrial biosphere model estimates. *J. Adv. Model. Earth Syst.* **12**, e2019MS001889 (2020).

60. Schewe, J. et al. State-of-the-art global models underestimate impacts from climate extremes. *Nat. Commun.* **10**, 1005 (2019).

61. Gentine, P. et al. Coupling between the terrestrial carbon and water cycles - A review. *Environ. Res. Lett.* **14**, 83003 (2019).

62. Bastos, A. et al. Impacts of extreme summers on European ecosystems: a comparative analysis of 2003, 2010 and 2018. *Philos. Trans. R. Soc. B Biol. Sci.* **375**, 20190507 (2020).

63. Fan, Y., Miguel-Macho, G., Jobbágy, E. G., Jackson, R. B. & Otero-Casal, C. Hydrologic regulation of plant rooting depth. *Proc. Natl Acad. Sci. USA* **114**, 10572–10577 (2017).

64. Haverd, V. et al. A new version of the CABLE land surface model (Subversion revision r4601) incorporating land use and land cover change, woody vegetation demography, and a novel optimisation-based approach to plant coordination of photosynthesis. *Geosci. Model Dev.* **11**, 2995–3026 (2018).

65. Melton, J. R. & Arora, V. K. Competition between plant functional types in the Canadian Terrestrial Ecosystem Model (CTEM) v. 2.0. *Geosci. Model Dev.* **9**, 323–361 (2016).

66. Lawrence, D. M. et al. The community land model version 5: description of new features, benchmarking, and impact of forcing uncertainty. *J. Adv. Model. Earth Syst.* **11**, 4245–4287 (2019).

67. Tian, H. et al. North American terrestrial CO₂ uptake largely offset by CH₄ and N₂O emissions: toward a full accounting of the greenhouse gas budget. *Clim. Change* **129**, 413–426 (2015).

68. Meiyappan, P., Jain, A. K. & House, J. I. Increased influence of nitrogen limitation on CO₂ emissions from future land use and land use change. *Glob. Biogechem. Cycles* **29**, 1524–1548 (2015).

69. Mauritsen, T. et al. Developments in the MPI-M Earth system model version 1.2 (MPI-ESM1.2) and its response to increasing CO₂. *J. Adv. Model. Earth Syst.* **11**, 998–1038 (2019).

70. Poulter, B., Frank, D. C., Hodson, E. L. & Zimmermann, N. E. Impacts of land cover and climate data selection on understanding terrestrial carbon dynamics and the CO₂ airborne fraction. *Biogeosciences* **8**, 2027–2036 (2011).

71. Lienert, S. & Joos, F. A Bayesian ensemble data assimilation to constrain model parameters and land-use carbon emissions. *Biogeosciences* **15**, 2909–2930 (2018).

72. Zaehle, S. & Friend, A. D. Carbon and nitrogen cycle dynamics in the O-CN land surface model: I. Model description, site-scale evaluation, and sensitivity to parameter estimates. *Glob. Biogechem. Cycles* **24**, 1–13 (2010).

73. Goll, D. S. et al. A representation of the phosphorus cycle for ORCHIDEE (revision 4520). *Geosci. Model Dev.* **10**, 3745–3770 (2017).
74. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Glob. Biogeochem. Cycles* **19**, 1–33 (2005).

75. Walker, A. P. et al. The impact of alternative trait-scaling hypotheses for the maximum photosynthetic carboxylation rate (Vcmax) on global gross primary production. *N. Phytol.* **215**, 1370–1386 (2017).

76. Joetzjer, E. et al. Improving the ISBACC land surface model simulation of water and carbon fluxes and stocks over the Amazon forest. *Geosci. Model Dev.* **8**, 1709–1727 (2015).

77. Kato, E., Kinoshita, T., Ito, A., Kawamiya, M. & Yamagata, Y. Evaluation of spatially explicit emission scenario of land-use change and biomass burning using a process-based biogeochemical model. *J. Land Use Sci.* **8**, 104–122 (2013).

78. Wei, Y. et al. The North American carbon program multi-scale synthesis and terrestrial model intercomparison project – Part 2: environmental driver data. *Geosci. Model Dev.* **7**, 2875–2893 (2014).

79. Dlugokencky, E. J., Thoning, K. W., Lang, P. M. & Tans, P. P. NOAA greenhouse gas reference from atmospheric carbon dioxide dry air mole fractions from the NOAA ESRL carbon cycle cooperative global air sampling network. ftp://aftp.cmdl.noaa.gov/data/trace_gases/co2/flask/surface/ (2017).

## Acknowledgements

We thank Dr. Martin Jung for his helpful comments on the manuscript and help in data preparation. We would like to thank the TRENDYv7 modelers for their simulations. We also thank Dr. Pieter Tans and Dr. Ed Dlugokencky for providing the CO₂ mole fraction data. S.P. acknowledges the support from the National Natural Science Foundation of China (41988101), the Second Tibetan Plateau Scientific Expedition and Research Program (Grant No. 2019QZKK0208), and the Xplorer Prize. A.B., P.C., F.C., M.O., and S.S. acknowledge the support from the European Space Agency Climate Change Initiative RECCAP2 project (contract no. 4000123002/18/I-NB). V.H. acknowledges the support from the Swiss National Science Foundation (grants no. P400P2_180784 and P4P4P2_194464).

## Author contributions

S.P., P.C., A.B., and X.W. designed the study. K.W. performed the analysis. K.W., S.P., A.B., and P.C. drafted the paper with contributions by X.W., C.R., P.G., F.C., V.W.H.,

C.H., M.O., S.I.S., and S.S. P.C. and A.B. contributed equally to this work. All authors contributed to the interpretation of the results.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-31175-w.

**Correspondence** and requests for materials should be addressed to Shilong Piao.

**Peer review information** *Nature Communications* thanks Kirsten Thonicke and the anonymous reviewer(s) for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022