# Compensatory water effects link yearly global land CO2 sink changes to temperature

Martin Jung, Markus Reichstein, Chris Schwalm, Chris Huntingford, Stephen Sitch, Anders Ahlström, Almut Arneth, Gustau Camps-Valls, Philippe Ciais, Pierre Friedlingstein, et al.

► To cite this version:

Martin Jung, Markus Reichstein, Chris Schwalm, Chris Huntingford, Stephen Sitch, et al. Compensatory water effects link yearly global land CO2 sink changes to temperature. Nature, 2017, 541 (7638), pp.516 - 520. 10.1038/nature20780. cea-01893465

HAL Id: cea-01893465

https://cea.hal.science/cea-01893465v1

Submitted on 14 Nov 2024

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>
**TITLE**

Compensatory water effects link yearly global land CO2 sink changes to temperature.

**AUTHORS**

Jung, M; Reichstein, M; Schwalm, CR; et al.

**JOURNAL**

Nature

**DEPOSITED IN ORE**

10 March 2017

This version available at

http://hdl.handle.net/10871/26422
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

# Compensatory water effects link yearly global land CO₂ sink changes to temperature

Martin Jung¹, Markus Reichstein¹,², Christopher R Schwalm³, Chris Huntingford⁴, Stephen Sitch⁵,

Anders Ahlström⁶,⁷, Almut Arneth⁸, Gustau Camps-Valls⁹, Philippe Ciais¹⁰, Pierre Friedlingstein¹¹,

Fabian Gans¹, Kazuhito Ichii¹²,¹³, Atul K. Jain¹⁴, Etsushi Kato¹⁵, Dario Papale¹⁶, Ben Poulter¹⁷, Botond

Raduly¹⁶,²⁰, Christian Rödenbeck¹⁸, Gianluca Tramontana¹⁶, Nicolas Viovy¹⁰, Ying-Ping Wang¹⁹, Ulrich

Weber¹, Sönke Zaehle¹,², Ning Zeng²¹

¹Max Planck Institute for Biogeochemistry, Department of Biogeochemical Integration, 07745 Jena, Germany

²Michael-Stifel-Center Jena for Data-driven and Simulation Science, 07743 Jena, Germany

³Woods Hole Research Center, 149 Woods Hole Rd, Falmouth, MA 02540, USA

⁴Centre for Ecology and Hydrology, Wallingford, Oxfordshire, OX10 8BB, U.K.

⁵College of Life and Environmental Sciences, University of Exeter, Exeter, UK.

⁶Department of Earth System Science, School of Earth, Energy and Environmental Sciences, Stanford University, Stanford, CA 94305, USA

⁷Lund University, Department of Physical Geography and Ecosystem Science, 223 62 Lund, Sweden.

⁸Karlsruhe Institute of Technology (KIT), Institute of Meteorology and Climate Research, 82467 Garmisch-Partenkirchen, Germany

⁹Image Processing Laboratory (IPL), C/ Catedrático José Beltran, 2. 46980 Paterna, València. Spain.

¹⁰Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ, 91191 Gif-sur-Yvette, France

¹¹College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QE, UK.

¹²Department of Environmental Geochemical Cycle Research, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173-25, Showa-machi, Kanazawa-ku, Yokohama, 236-0001, Japan

¹³Center for Global Environmental Research, National Institute for Environmental Studies, 16-2, Onogawa, Tsukuba, 305-8506, Japan
31 <sup>14</sup>Department of Atmospheric Sciences, University of Illinois, Urbana, IL 61801, USA

32 <sup>15</sup>Global Environment Programs, The Institute of Applied Energy (IAE), Tokyo, 105-0003, Japan

33 <sup>16</sup>Department for Innovation in Biological, Agro-food and Forest systems (DIBAF), Via San Camillo de Iellis snc, 01100 Viterbo, Italy.
34

35 <sup>17</sup>Department of Ecology, Montana State University, Bozeman, MT 59717

36 <sup>18</sup>Max Planck Institute for Biogeochemistry, Department of Biogeochemical Systems, 07745 Jena,
37 Germany

38 <sup>19</sup>CSIRO Ocean and Atmosphere, PMB #1, Aspendale, Victoria 3195, Australia

39 <sup>20</sup>Dept. of Bioengineering, Sapientia Hungarian University of Transylvania, 530104 M-Ciuc, Romania

40 <sup>21</sup>Dept. of Atmospheric and Oceanic Science, University of Maryland, College Park, MD 20742, USA

41
Large interannual variations in the measured growth rate of atmospheric carbon dioxide originate primarily from fluctuations in the carbon uptake by land ecosystems¹⁻³. It remains uncertain, however, to what extent temperature and water availability control the carbon balance of land ecosystems across spatial and temporal scales³⁻¹⁴. Here we use eddy covariance data-derived empirical models¹⁵ and process based models¹⁶,¹⁷ to investigate the effect of changes in temperature and water availability on gross primary productivity (GPP), terrestrial ecosystem respiration (TER) and net ecosystem exchange (NEE) at local and global scales. We find that water availability is the predominant driver of the interannual variability in GPP, TER and, to a lesser extent, NEE at the local scale. When integrated globally, however, temporal NEE variability is mostly driven by temperature fluctuations (R²≥0.84). We suggest that this apparent paradox can be explained by two compensatory water effects. Temporal water driven GPP and TER variations compensate locally, dampening water-driven NEE variability. Spatial water availability anomalies also compensate, leaving a dominant temperature signal in the year-to-year fluctuations of the land carbon sink. These findings help reconcile seemingly contradictory reports regarding the importance of temperature and water in controlling the interannual variability of the terrestrial carbon balance³⁻⁶,⁹,¹¹,¹²,¹⁴. Our study indicates that spatial climate co-variation drives the global carbon cycle response.

Large interannual variations in the recent measured atmospheric CO₂ growth rates originate primarily from fluctuations in carbon uptake by land ecosystems, rather than from oceans or variations in anthropogenic emissions¹⁻³. There is a general consensus that the tropical region contributes the most to terrestrial carbon variability¹,⁸,¹⁸,¹⁹. The observed positive correlation between mean tropical land temperature and CO₂ growth rate³,⁵,⁶,¹²,¹³ implies smaller land carbon uptake and enhanced atmospheric CO₂ growth during warmer years with a sensitivity of about 5 GtC yr⁻¹K⁻¹. There is a tight relationship between this sensitivity on interannual time scales and long-term changes in terrestrial carbon per degree of warming across multiple climate carbon-cycle models⁶. Despite this strong emergent relationship with mean tropical land temperature, several studies suggest that variations in water availability play an important⁸,¹⁰,¹¹,¹⁴, even a dominant role⁴,⁹, in shaping the interannual variability of the carbon balance of extensive semi-arid and sub-tropical systems. Furthermore, the recent doubling of the tropical carbon cycle sensitivity to interannual temperature variability has been linked to interactions with changing moisture regimes¹³. A full understanding of the processes governing the climatic controls of terrestrial carbon cycling on interannual time scales and across spatial scales is therefore still lacking. Here we show that the “temperature vs. water” debate can be resolved by simultaneously assessing the carbon cycle response to fluctuations in both temperature and water availability at both local and global scales.

Using both machine learning algorithms and process-based global land models, we derived spatial and temporal patterns of the interannual variability (IAV) of CO₂ uptake by plants via photosynthesis (gross primary production, GPP) and of CO₂ loss through respiration (terrestrial ecosystem respiration, TER). This allows analysis of net CO₂ ecosystem exchange (NEE=TER-GPP) IAV. Machine learning algorithms were used to translate gridded inputs of daily air temperature, water availability and radiation, among others¹⁵, into time varying 0.5° grids of TER and GPP for the 1980-2013 period (FLUXCOM, see Methods). Three machine learning algorithms were trained on FLUXNET²⁰ based in situ TER and GPP flux estimates from two flux partitioning methods²¹,²². These three fitting algorithms combined with two partitioning methods provided six sets of GPP and TER estimates each, which combined yield 36 FLUXCOM NEE ensemble members. In a complementary approach,
we examined simulations of GPP and TER from an ensemble of seven global land surface or dynamic vegetation models¹⁶,¹⁷ (TRENDYv3, see Methods). These process-based model simulations follow a common protocol and used the same climate forcing data set as the observation-based FLUXCOM models. Both sets of results are expected to be more uncertain in the tropics due to less reliable climate and satellite based inputs and a sparse coverage of flux measurements²³.

We analysed FLUXCOM and TRENDYv3 simulations independently, but in a consistent manner. We derived NEE as the difference between TER and GPP, i.e., a positive value of NEE indicates a flux of carbon from the land to the atmosphere. To isolate IAV we detrended GPP and TER for each grid cell and month (see Methods). We find that global patterns of NEE interannual variability are consistent between FLUXCOM and TRENDYv3 (EDF 1, SI-1). Both approaches reproduce (r ~ 0.8) the globally integrated NEE IAV derived from atmospheric CO₂ concentration measurements and transport²⁴. Both approaches also show the largest IAV in the tropics (EDF 1). To obtain the contributions of different environmental variables to IAV, we decomposed carbon flux anomalies (ΔFlux) of each year (y), month (m), and grid cell (s) into their additive components forced by detrended anomalies of temperature (ΔTEMP), shortwave incoming radiation (ΔRAD), and soil-moisture related water availability (ΔWAI, see Methods):

$$ 
\begin{array}{l}
\Delta Flux_{s,m,y} = a_{s,m}^{TEMP} \times \Delta TEMP_{s,m,y} + a_{s,m}^{RAD} \times \Delta RAD_{s,m,y} + a_{s,m}^{WAI} \times \Delta WAI_{s,m,y} + \varepsilon_{s,m,y} \\
\Delta Flux_{s,m,y} \approx \Delta Flux_{s,m,y}^{TEMP} + \Delta Flux_{s,m,y}^{RAD} + \Delta Flux_{s,m,y}^{WAI}.
\end{array}
\qquad \text{EQ (1)}
$$

Here $a_{s,m}$ represents the estimated sensitivity of the flux anomaly, $\Delta Flux_{s,m,y}$ (GPP or TER) to each respective climate forcing anomaly ($\Delta TEMP$, $\Delta RAD$, $\Delta WAI$) for a given grid cell and month, and $\varepsilon_{s,m,y}$ is the residual error term. The product of a given sensitivity (e.g. $a^{TEMP}$) and corresponding climate forcing anomaly (e.g. $\Delta TEMP$) constitutes the flux anomaly component driven by this climate factor (e.g. GPP$^{TEMP}$). Thus, Eq.1 estimates the contributions of temperature, radiation, and water availability anomalies to carbon flux anomalies (see SI-2 for verification).

Our analysis reveals a contrasting pattern of NEE IAV controlled by temperature or moisture, depending on spatial scale. At the global scale, temperature drives spatially-integrated NEE IAV (Fig.1 a,b, compare green and black curves), in line with previous findings based on correlations between anomalies in temperature and CO₂ growth rate³,⁵,⁶,¹²,¹³. Globally integrated NEE anomalies due to variations in radiation (NEE$^{WAI}$) and water availability (NEE$^{WAI}$) play only a minor role (compare blue and black curves in Fig. 1a,b). The dominant global influence of temperature is in contrast to the dominant local influence of water availability when analyzing all grid cells individually (Fig 1 c,d, zonal mean of grid cell IAV; compare blue and black curves). Radiation causes the smallest NEE IAV at grid cell level (red curve in Fig.1c,d) but there are indications based on other climate forcing data that radiation could play a more important role than temperature locally (SI-3). Temperature variations are important for NEE IAV (green curve in Fig.1c,d) in high latitudes and the inner tropics, but in general, the grid cell average water related NEE variability (NEE$^{WAI}$, blue curve) is larger. Water related NEE variability peaks at subtropical latitudes where semi-arid ecosystems dominate. This finding is consistent with studies emphasizing the role of water limited semi-arid ecosystems on global NEE IAV⁴,⁹. We now assess how this can be reconciled with the emergent temperature control of globally integrated NEE IAV. Going from grid-cell to global scale shifts the emerging controls on NEE IAV from water availability (local) towards temperature (global).
127

[Insert Figure 1 around here]

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

We hypothesized that the dominance of temperature in globally integrated NEE IAV results from a stronger compensation of positive and negative NEE<sup>WAI</sup> anomalies between different grid cells compared to NEE<sup>TEMP</sup> when going from local to global scale. To test this, we first illustrate the dominant spatial patterns of temperature vs. water compensation using empirical orthogonal functions (EOF) of the annual NEE<sup>TEMP</sup> and NEE<sup>WAI</sup> anomalies (Fig. 2 a-d). Here, the leading EOF of NEE<sup>WAI</sup> (~10% variance explained) has strong anti-correlated spatial patterns of positive and negative values (Fig 2c,d), which correspond to ENSO imprints on moisture effects (R<sup>2</sup> with Nino 3.4 SST index<sup>25</sup> of 0.75). In comparison, the leading EOF of NEE<sup>TEMP</sup> (~22% variance explained) shows a more spatially uniform response, and in particular across the tropics (Fig 2a,b). This pattern of much larger spatial coherence of NEE<sup>TEMP</sup> anomalies, compared to NEE<sup>WAI</sup> anomalies, is also evident in their respective sums of positive and negative covariances among all grid cells (inset pie charts in Fig 2. a-d). For NEE<sup>TEMP</sup> the sum of positive covariances is far larger than the negative ones (79% vs. 21%), whereas positive and negative covariances are almost in balance (53% vs. 47%) for NEE<sup>WAI</sup>. As a consequence of the larger spatial coherence of NEE<sup>TEMP</sup> anomalies, as compared to NEE<sup>WAI</sup> anomalies, we observe a shift of the dominant NEE IAV control from water at the local scale to temperature at the global scale. We illustrate this change in Fig 2e,f by presenting relative dominance of water and temperature related NEE IAV for increasing levels of spatial aggregation. This is a robust feature within and among FLUXCOM and TRENDY approaches (EDF 2). We also find that the rise and decay of NEE<sup>TEMP</sup> and NEE<sup>WAI</sup> dominance respectively with spatial scale occurs in all major biomes (SI-4). This pattern is likely related to the different climatic characteristics of precipitation and air temperatures, with the former, but not the latter, being associated with moisture conservation and offsetting spatial anomaly patterns.

150

[Insert Figure 2 around here]

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

We now proceed to assess how local water and temperature related NEE IAV emerges from the interaction of photosynthesis (GPP) and respiration (TER) processes. We compare the magnitudes of water vs. temperature driven GPP and TER variability and find that WAI is overall the most important factor controlling local IAV of both gross fluxes (Fig. 3 a-d), with particularly large variability in both fluxes in semi-arid regions (SI-4, 5). However, the local IAV of NEE related to WAI (NEE<sup>WAI</sup>, Fig. 3e, f) is reduced compared to the components GPP<sup>WAI</sup> and TER<sup>WAI</sup>. Our results indicate that, in addition to the spatial compensation of NEE<sup>WAI</sup>, discussed above, there is also a local compensation mechanism, whereby GPP<sup>WAI</sup> and TER<sup>WAI</sup> co-vary and thus locally counterbalance each other (Fig. 4 a, b). This is likely due to the concomitant positive relationship of soil moisture with productivity and with respiration. The combined effect is a smaller net effect of WAI on NEE. Specifically, two thirds of the WAI effect on GPP is offset by the WAI effect on TER (0.67±0.33 for FLUXCOM, 0.69±0.14 for TRENDY; mean slope ± s.d. across ensemble members of global TER<sup>WAI</sup> vs. GPP<sup>WAI</sup>). These patterns are qualitatively consistent between the data-driven FLUXCOM (Fig. 4) and process-based TRENDY models (EDF. 3) and agree with previous observations of simultaneous declines of GPP and TER<sup>26-3025-29</sup> during droughts. However, magnitudes of TER<sup>WAI</sup> vs. GPP<sup>WAI</sup> covariances differ substantially among model ensemble members (EDF 4). This likely reflects large uncertainty of respiration processes to moisture variations while flux partitioning uncertainties seem negligible (SI-6).

168

[Insert Figure 3 around here]
169 In contrast to offsetting NEE water effects, our analysis indicates a weak local temperature amplification effect of GPP and TER IAV in the tropics. Local temperature effects on GPP and TER IAV are inversely correlated over the tropics (Fig. 4d). This is because GPP decreases with increasing temperature, likely due to the exceedance of the thermal optimum of photosynthesis, whereas respiration increases with temperature. Thus increasing temperatures in the tropics reduce NEE by reducing GPP and increasing TER. However, due to lower variances of the temperature components of GPP and TER (Fig. 3a-d), this local temperature amplification effect in the tropics is quantitatively negligible (Fig. 4c) compared to the local water compensation effect (Fig. 4d). Overall, this causes the difference of temperature vs. water forced variability of NEE to be smaller compared to the influence of these drivers on the gross fluxes (compare distance between blue and green curves in Fig. 3 a-d vs. e, f).

180

[Insert Figure 4 around here]

181 Our analysis shows water availability as the overall dominant driver of the interannual variability of photosynthesis and respiration at local scales, even though this water signal is effectively absent in the globally integrated NEE interannual variability. This pattern is driven by: 1) the local compensatory effects of water availability on GPP and TER, and 2) the spatial anti-correlation of water controlled NEE anomalies; and thus a compensation in space. These two compensatory water effects leave temperature as the dominant factor globally, which resolves why there have been conflicting conclusions surrounding whether NEE interannual variability is forced thermally or hydrologically. Our analysis implies that climate does not only force the carbon cycle locally, but that, perhaps more importantly, the spatial covariation of climate variables drives the integrated global carbon cycle response. Consequently, any analysis conducted on integrated signals over larger regions precludes inferences on the driving mechanisms at ecosystem scale. Likewise, the apparent temperature dominated interannual variability of the residual land sink, a traditional target of global carbon cycle modelers, contains little information on local carbon cycle processes. Our findings suggest that potential changes in spatial covariations among climate variables associated with global change may drive apparent changes of carbon cycle sensitivities and perhaps even the strength of climate-carbon cycle feedbacks.

197

198 References

199

200 1 Bousquet, P. et al. Regional changes in carbon dioxide fluxes of land and oceans since 1980. *Science* **290**, 1342-1346 (2000).

201 2 Keeling, C. D., Piper, S. C. & Heimann, M. in *Aspects of Climate Variability in the Pacific and Western Americas* (ed D. H. Peterson) 305-363 (Washington, D.C., 1989).

202 3 Kindermann, J., Wurth, G., Kohlmaier, G. H. & Badeck, F. W. Interannual variation of carbon exchange fluxes in terrestrial ecosystems. *Global Biogeochemical Cycles* **10**, 737-755, doi:10.1029/96gb02349 (1996).

207 4 Ahlstrom, A. et al. The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink. *Science* **348**, 895-899, doi:10.1126/science.aaa1668 (2015).

209 5 Anderegg, W. R. L. et al. Tropical nighttime warming as a dominant driver of variability in the terrestrial carbon sink. *Proceedings of the National Academy of Sciences* **112**, 15591-15596, doi:10.1073/pnas.1521479112 (2015).
<table>
  <tr>
    <td>212</td>
    <td>6</td>
    <td>Cox, P. M. et al. Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability. <i>Nature</i> 494, 341-344, doi:10.1038/nature11882 (2013).</td>
  </tr>
  <tr>
    <td>213</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>214</td>
    <td>7</td>
    <td>Nemani, R. R. et al. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. <i>Science</i> 300, 1560-1563 (2003).</td>
  </tr>
  <tr>
    <td>216</td>
    <td>8</td>
    <td>Peylin, P. et al. Multiple constraints on regional CO2 flux variations over land and oceans. <i>Global Biogeochemical Cycles</i> 19, doi:10.1029/2003gb002214 (2005).</td>
  </tr>
  <tr>
    <td>217</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>218</td>
    <td>9</td>
    <td>Poulter, B. et al. Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle. <i>Nature</i> 509, 600-603 (2014).</td>
  </tr>
  <tr>
    <td>220</td>
    <td>10</td>
    <td>Tian, H. Q. et al. Effect of interannual climate variability on carbon storage in Amazonian ecosystems. <i>Nature</i> 396, 664-667, doi:10.1038/25328 (1998).</td>
  </tr>
  <tr>
    <td>221</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>222</td>
    <td>11</td>
    <td>Wang, J., Zeng, N. & Wang, M. Interannual variability of the atmospheric CO2 growth rate: roles of precipitation and temperature. <i>Biogeosciences</i> 13, 2339-2352, doi:10.5194/bg-13-2339-2016 (2016).</td>
  </tr>
  <tr>
    <td>225</td>
    <td>12</td>
    <td>Wang, W. et al. Variations in atmospheric CO2 growth rates coupled with tropical temperature. <i>Proceedings of the National Academy of Sciences of the United States of America</i> 110, 13061-13066, doi:10.1073/pnas.1219683110 (2013).</td>
  </tr>
  <tr>
    <td>226</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>227</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>228</td>
    <td>13</td>
    <td>Wang, X. et al. A two-fold increase of carbon cycle sensitivity to tropical temperature variations. <i>Nature</i> 506, 212-215 (2014).</td>
  </tr>
  <tr>
    <td>229</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>230</td>
    <td>14</td>
    <td>Zeng, N., Mariotti, A. & Wetzel, P. Terrestrial mechanisms of interannual CO(2) variability. <i>Global Biogeochemical Cycles</i> 19, doi:10.1029/2004gb0022763 (2005).</td>
  </tr>
  <tr>
    <td>231</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>232</td>
    <td>15</td>
    <td>Tramontana, G. et al. Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms. <i>Biogeosciences</i> 13, 4291-4313, doi:10.5194/bg-13-4291-2016 (2016).</td>
  </tr>
  <tr>
    <td>233</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>234</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>235</td>
    <td>16</td>
    <td>Le Quéré, C. et al. Global carbon budget 2014. <i>Earth System Science Data Discussions</i> 7, 521-610 (2014).</td>
  </tr>
  <tr>
    <td>236</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>237</td>
    <td>17</td>
    <td>Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon dioxide. <i>Biogeosciences</i> 12, 653-679, doi:10.5194/bg-12-653-2015 (2015).</td>
  </tr>
  <tr>
    <td>238</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>239</td>
    <td>18</td>
    <td>Gurney, K. R., Baker, D., Rayner, P. & Denning, S. Interannual variations in continental-scale net carbon exchange and sensitivity to observing networks estimated from atmospheric CO(2) inversions for the period 1980 to 2005. <i>Global Biogeochemical Cycles</i> 22, doi:Gb302510.1029/2007gb003082 (2008).</td>
  </tr>
  <tr>
    <td>240</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>241</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>242</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>243</td>
    <td>19</td>
    <td>Peylin, P. et al. Global atmospheric carbon budget: results from an ensemble of atmospheric CO2 inversions. <i>Biogeosciences</i> 10, 6699-6720, doi:10.5194/bg-10-6699-2013 (2013).</td>
  </tr>
  <tr>
    <td>244</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>245</td>
    <td>20</td>
    <td>Baldocchi, D. TURNER REVIEW No 15. 'Breathing' of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. <i>Australian Journal of Botany</i> 56, 1-26 (2008).</td>
  </tr>
  <tr>
    <td>246</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>247</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>248</td>
    <td>21</td>
    <td>Lasslop, G. et al. Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation. <i>Global Change Biology</i> 16, 187-208 (2010).</td>
  </tr>
  <tr>
    <td>249</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>250</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>251</td>
    <td>22</td>
    <td>Reichstein, M. et al. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. <i>Global Change Biology</i> 11, 1424-1439 (2005).</td>
  </tr>
  <tr>
    <td>252</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>253</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>254</td>
    <td>23</td>
    <td>Schimel, D. et al. Observing terrestrial ecosystems and the carbon cycle from space. <i>Global Change Biology</i> 21, 1762-1776, doi:10.1111/gcb.12822 (2015).</td>
  </tr>
  <tr>
    <td>255</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>256</td>
    <td>24</td>
    <td>Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO2 flux history 1982-2001 inferred from atmospheric data using a global inversion of atmospheric transport. <i>Atmos. Chem. Phys.</i> 3, 1919-1964, doi:10.5194/acp-3-1919-2003 (2003).</td>
  </tr>
  <tr>
    <td>257</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>258</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>259</td>
    <td>25</td>
    <td>Rayner, N. A. et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. <i>Journal of Geophysical Research-Atmospheres</i> 108, doi:10.1029/2002jd002670 (2003).</td>
  </tr>
  <tr>
    <td>260</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>261</td>
    <td></td>
    <td></td>
  </tr>
</table>
262 26 Biederman, J. A. et al. Terrestrial carbon balance in a drier world: the effects of water availability in southwestern North America. *Global Change Biology* **22**, 1867-1879, doi:10.1111/gcb.13222 (2016).

265 27 Bonal, D., Burban, B., Stahl, C., Wagner, F. & Hérault, B. The response of tropical rainforests to drought—lessons from recent research and future prospects. *Annals of Forest Science* **73**, 27-44, doi:10.1007/s13595-015-0522-5 (2016).

268 28 Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat and drought in 2003. *Nature* **437**, 529-533 (2005).

270 29 Schwalm, C. R. et al. Assimilation exceeds respiration sensitivity to drought: A FLUXNET synthesis. *Global Change Biology* **16**, 657-670 (2010).

272 30 van der Molen, M. K. et al. Drought and ecosystem carbon cycling. *Agricultural and Forest Meteorology* **151**, 765-773 (2011).

274

## Acknowledgements

276 We thank Philippe Peylin for providing RECCAP inversion results. We acknowledge Paul Bodesheim for help with the mathematical notations, Jacob Nelson for proof reading the SI, Silvana Schott for help with art work, and Gerhard Boenisch, Linda Maack, and Peer Koch for help on archiving the FLUXCOM data. MJ, MR, DP acknowledge funding from the EU FP7 project GEOCARBON (grant no. 283080) and the EU H2020 BACI project (grant no. 640176). FG, and MR acknowledge the European Space Agency ESA for funding the "Coupled Biosphere-Atmosphere virtual LABoratory, CAB-LAB'. SZ acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (QUINCY; grant no. 647204). AAr acknowledges support from the EU FP7 project LUC4C (grant no. 603542). CRS was supported by National Aeronautics and Space Administration (NASA) grants NNX12AK12G, NNX12AP74G, NNX10AG01A, and NNX11AO08A. PC acknowledges support from the European Research Council Synergy grant ERC-2013-SyG-610028 IMBALANCE-P. SS acknowledges the support of the Natural Environment Research Council (NERC) South AMerican Biomass Burning Analysis (SAMBBA) project grant code NE/J010057/1. CH is grateful for support from the NERC CEH National Capability fund. AAh acknowledge support from The Royal Physiographic Society in Lund (Birgit and Hellmuth Hertz' Foundation) and the Swedish Research Council (637-2014-6895). GCV was supported by the EU under the European Research Council (ERC) consolidator grant SEDAL-647423

293 Author Contributions

294 MJ and MR designed the analysis. MJ carried out the analysis and wrote the manuscript with contributions from all authors. MJ, CRS, GCV, FG, KI, DP, BR, GT, and UW contributed to FLUXCOM results. SS, PF, CH, AAI, Aar, PC, AKJ, EK, BP, NV, YPW, and NZ contributed to TRENDY results.

297 Author Information

298 Reprints and permissions information is available at http://www.nature.com/reprints. The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to mjung@bgc-jena.mpg.de.

301

302
303 **Figure captions**

304

305 **Figure 1: Climatic controls on NEE IAV at global and local scales for the period 1980-2013 derived from machine learning based (FLUXCOM) and process-based (TRENDY) models.** *The comparison of globally integrated annual NEE anomalies with NEE anomalies driven only by temperature, water availability, and radiation (a, b) shows temperature as dominant global control. R² values between the climatic NEE components and total NEE are given in the respective colour. Mean grid cell IAV magnitude (see Equation 3 in Methods) in panels (c) and (d) of NEE components for latitudinal bands shows water as dominant local control. Uncertainty bounds where given as shaded area reflect the spread among FLUXCOM or TRENDY ensemble members (±1 s.d.).*

313 **Figure 2: Effects of spatial co-variation and scale on temperature vs. water control of NEE IAV for FLUXCOM and TRENDY models.** *Spatial patterns of the first empirical orthogonal function of annual NEE<sup>TEMP</sup> (a, b), and NEE<sup>WAI</sup> (c, d) anomalies (see Methods) show large spatial coherence for NEE<sup>TEMP</sup> (dominant positive values) and anti-correlated patterns for NEE<sup>WAI</sup> (positive and negative values; magnitudes are not informative and were omitted for clarity). This is underpinned in the inset pie charts which show the proportion of total positive (black) and negative (gray) co-variances among grid cells for NEE<sup>TEMP</sup> and NEE<sup>WAI</sup> anomalies (see Equation 4 and 5 in Methods). Panels e, f present how the relative dominance (see Equation 6 in Methods) of NEE<sup>TEMP</sup> (green) increases with successive spatial aggregation, while the relative dominance of NEE<sup>WAI</sup> (blue) decreases. Outer uncertainty bounds in e,f, given as shaded area refer to the spread among respective ensemble members (±1 s.d.); inner uncertainty bounds refer to ±1 s.d. with respect to the change of relative dominance with spatial aggregation (see Equation 7 in Methods).*

325 **Figure 3: Latitudinal patterns of water and temperature driven IAV of gross carbon fluxes (GPP and TER) and NEE for FLUXCOM and TRENDY models.** *IAV magnitude (see Equation 3 in Methods) of the WAI component is much larger than the IAV of the TEMP component for gross fluxes (a-d), while this difference is smaller for NEE due to compensation. Uncertainty bounds as shaded area reflect the spread among FLUXCOM or TRENDY ensemble members (±1 s.d.).*

330 **Figure 4: Spatial patterns of covariance and correlation of WAI and TEMP driven GPP and TER IAV for FLUXCOM models.** *Maps of the covariance of annual anomalies (see Equation 8 in Methods) of GPP and TER climatic components show large compensation effects (positive covariance) for WAI (a) but nearly no covariance for TEMP (c). Correlations between GPP<sup>WAI</sup> and TER<sup>WAI</sup> are large and ubiquitous positive (b) while correlations among GPP<sup>TEMP</sup> and TER<sup>TEMP</sup> are weaker with a distinct spatial pattern of negative correlations in hot regions (d). All results refer to the mean of all FLUXCOM ensemble members. See EDF 3 for equivalent TRENDY results, and EDF 4 for uncertainties.*

337
338 Methods

339 Global carbon flux data sets

340 FLUXCOM. Three machine learning methods were trained on daily carbon flux estimates from 224
341 flux tower sites using meteorological measurements and satellite data as inputs¹⁵: Random Forests³¹,
342 Artificial Neural Networks³², Multivariate Adaptive Regression Splines³³. Models were trained
343 separately for two variants of GPP and TER, derived from the flux partitioning methods of Reichstein
344 et al.²² and Lasslop et al.²¹. Each method used the same 11 input driver data listed in Table SI-7. This
345 set of driver data was obtained from an extensive variable selection analysis¹⁵,³⁴. Details along with
346 extensive model evaluation based on cross-validation are given in Tramontana et al.¹⁵.

347 To produce spatio-temporal grids of carbon fluxes, the trained machine learning algorithms require
348 only spatio-temporal grids of its input driver data³⁵. We forced the models with grids of 0.5° spatial
349 resolution and daily time step for the period 1980-2013³⁶. High-resolution satellite based predictor
350 variables (see Table SI-7) were tiled by plant functional type (PFT), i.e. grids for each PFT containing
351 the mean value per PFT and time step at 0.5° were created. The PFT distribution originates from the
352 majority class of annually resolved MODIS land cover product (collection 5)³⁷ for each high-resolution
353 pixel. Climatic predictor variables are based on CRUNCEPv6
354 (http://esgf-extra.cea.fr/thredds/catalog/store/p529viov/cruncep/V6_1901_2014/catalog.html) to
355 be consistent with the TRENDY ensemble. CRUNCEPv6 is based on a merged product of Climate
356 Research Unit (CRU) observation based monthly 0.5° climate variables³⁸ (1901 – 2013) and the high
357 temporal (6-hourly) resolution NCEP reanalysis. The variables affected by the climate forcing data set
358 are marked in Table SI-7. Among the 11 predictor variables, only temperature, radiation, and water
359 availability can generate interannual variability. The water availability index (WAI, see supplement 3
360 in Tramontana et al. ¹⁵) is based on a simple dynamic soil water balance model, which was driven
361 with daily precipitation and potential evapotranspiration by CRUNCEPv6 (see SI-8 for cross-
362 consistency with TRENDY based soil moisture). The machine learning models were run at for each
363 plant functional type (PFT) separately, and a weighted mean over the PFT fractions was obtained for
364 each grid-cell. The PFT distribution is representative of the period 2001-2012; no land cover change
365 was considered. Empirical models were run to spatially estimate GPP and TER. Then NEE was derived
366 by the carbon mass balance approach (NEE = TER-GPP), which allows for decomposing precisely of
367 how NEE IAV emerges from (co-)variations of TER and GPP. We verify that NEE IAV derived as TER-
368 GPP is consistent with upscaling NEE directly (SI-6). Overall 36 combinations of NEE were derived by
369 considering all possible combinations of TER-GPP realizations resulting from different machine
370 learning approaches, and flux partitioning variants. The individual model runs were finally aggregated
371 to monthly means.

372 TRENDY. We used simulations of seven Dynamic Global Vegetation Models (DGVMs) from the
373 TRENDY v3 ensemble¹⁶,¹⁷ for the period 1980-2013, which have a spatial resolution of 0.5° (model
374 simulations with coarser resolution were omitted): CABLE³⁹, ISAM⁴⁰, LPJ⁴¹, LPJ-GUESS⁴², ORCHIDEE⁴³,
375 VEGAS¹⁴, VISIT⁴⁴. These models were forced by a common set of input datasets and experimental
376 protocol (experiment 'S2')¹⁶,¹⁷. Climate forcing (CRUNCEPv6) is the same as for FLUXCOM. Global
377 atmospheric CO₂ was derived from ice core and NOAA monitoring station data, and provided at
378 annual resolution over the period 1860-2013¹⁶. DGVMs were run from preindustrial steady state
379 (NEE = 0) with changing fields of climate and atmospheric CO₂ concentration over the 20thC. Land
380 Use and Land cover changes were not considered. For consistency with FLUXCOM, NEE was derived
as the difference between terrestrial ecosystem respiration (TER) and GPP, i.e. fire emissions available from some models were not included. Terrestrial ecosystem respiration was calculated as the sum of simulated autotrophic and heterotrophic respiration.

384

385 **Analysis**

386 **Anomalies and decomposition.** Detrended monthly anomalies were obtained by removing the linear trend over years for each pixel and month (least squares fitting), which also centers the mean to zero 387 for a given pixel and month. This procedure was applied consistently to GPP, and TER, shortwave 388 radiation (RAD), air temperature (TEMP), and water availability (WAI), FLUXCOM and TRENDY 389 simulations. For TRENDY models the simulated soil moisture was used instead of WAI. The resulting 390 IAV of GPP and TER was decomposed into the contributions forced by TEMP, RAD, and WAI following 391 Eq.1 using a multiple linear (ordinary least squares) regression with zero intercept for each pixel and 392 month. NEE sensitivities and NEE components were derived from GPP and TER results, which is 393 equivalent to decomposing NEE (=TER-GPP) directly. We validate and discuss the approximation of 394 IAV contributions by Eq.1 in SI-2. 395

396 **Notations.** All analysis is based on detrended monthly anomalies (Eq. 1) aggregated to annual means. 397 For simplicity, we omit the Δ notation for ‘anomaly’ in the following. Superscripts ‘TEMP’, ‘WAI’, 398 ‘RAD’ refer to surface air temperature, water availability, and incoming shortwave radiation of a 399 respective carbon flux anomaly. Subscripts ‘s’, ‘y’, ‘e’ refer to indexes of grid cell, year, and ensemble 400 member respectively. The mean and standard deviation are denoted as μ and σ respectively, where 401 the subscripts of these operators tell whether the operation is done over grid cells (e.g. μs is an 402 average over all grid cells), years (e.g. σy is the standard deviation over the years), or ensemble 403 members. All main results refer to the mean of FLUXCOM or TRENDY ensemble members (μe) and 404 the standard deviation (σe) is used as uncertainty estimate. Whenever we calculated a mean over 405 0.5° grid cells (μs) we accounted for different grid cell areas (area weighted mean) and used a 406 consistent mask of valid values between FLUXCOM and TRENDY. Because several analyses are 407 referenced with respect to the sum of climatic components of NEE we denote NEE*: 408

$$
408 \qquad NEE_{s,y}^* = NEE_{s,y}^{TEMP} + NEE_{s,y}^{WAI} + NEE_{s,y}^{RAD} \qquad \text{EQ (2)}
$$

409 **Spatial patterns of IAV magnitude (e.g. Fig. 1c,d & 3).** To describe spatial patterns of IAV magnitude 410 (M) of climatic components of carbon fluxes (e.g. GPP<sup>WAI</sup>) we computed the standard deviation of its 411 annual values (σ<sub>y</sub>) for each grid cell (s). This standard deviation was then normalized by the mean (μ<sub>s</sub>) 412 temporal standard deviation (σ<sub>y</sub>) of NEE* to provide a relative metric of IAV magnitude, where values 413 above 1 indicate IAV magnitudes larger than average NEE* IAV. This scaling accounts for the known 414 underestimation of IAV magnitude in the upscaling approach<sup>35</sup> but does not change any patterns. 415

415

$$
416 \qquad M_s = \frac{\sigma_y(\text{Flux}_{s,y}^{\text{COMP}})}{\mu_s(\sigma_y(NEE_{s,y}^*))} \qquad \text{EQ (3)}
$$

417 Fig. 1c,d shows mean and standard deviations across ensemble members (μe and σe) for NEE
418 components for latitudinal bins of 5°. The same holds for Fig.3 which shows also GPP and TER
419 components.
420 **Empirical orthogonal functions and spatial covariances (Fig. 2a-d).** We first calculated mean spatio-
421 temporal grids of NEE climatic components across ensemble members ($\mu_e(NEE_{s,y,e}^{COMP})$). We then
422 multiplied those with grid cell areas to convert flux densities into fluxes per grid cell, and normalized
423 them by the standard deviation of NEE* across time and space ($\sigma_{s,y}(\mu_e(NEE_{s,y,e}^{*}))$)). Empirical
424 orthogonal functions were then computed for each climatic component without additional scaling in
425 MATLAB using the 'pca' function. The spatial pattern of first principle components (leading EOFs) of
426 NEE<sup>TEMP</sup> and NEE<sup>WAI</sup> was plotted with the same color scale. The values on the color bar themselves
427 are not informative and were therefore omitted for clarity. The leading EOF explains about 22% of
428 spatial NEE<sup>TEMP</sup> variance and ~10% of spatial NEE<sup>WAI</sup> variance in both FLUXCOM and TRENDY
429 ensemble means.

430 To quantify the degree of spatial covariance of NEE climatic components (inset pie charts in Fig. 2a-d)
431 we calculated a large covariance matrix of all grid cells vs all grid cells for each NEE climatic
432 component (annual anomalies multiplied with grid cell area), where the elements of this covariance
433 matrix ($c_{i,j}^{\text{COMP}}$) were calculated according to Equation (4):

$$
434 \qquad c_{i,j}^{\text{COMP}} = \text{cov}_y(NEE_{si,y}^{\text{COMP}}, NEE_{sj,y}^{\text{COMP}}) \qquad \text{EQ (4)}
$$

435 Here *i* and *j* index the two grid cells for which the covariance is calculated. By definition the variance
436 of the globally integrated anomalies equals the sum of all terms in the covariance matrix. To
437 determine the share of positive vs negative spatial covariance of the total variance, we summed
438 positive and negative covariance terms respectively (Equation 5). The sum of variances (the diagonal
439 of the covariance matrix where *i*=*j*) was omitted in the pie charts because they accounted for less
440 than 1% of the covariance budget.

$$
441 \quad t\operatorname{cov}_+^{\mathrm{COMP}} = \sum_{i=1} \sum_{j \neq i} c_{i,j}^{\mathrm{COMP}} |c_{i,j}^{\mathrm{COMP}} > 0; t\operatorname{cov}_-^{\mathrm{COMP}} = \sum_{i=1} \sum_{j \neq i} c_{i,j}^{\mathrm{COMP}} |c_{i,j}^{\mathrm{COMP}} < 0 \quad \text{EQ (5)}
$$

442 **Scale dependence of relative dominance of NEE<sup>TEMP</sup> and NEE<sup>WAI</sup> (Fig. 2e,f).** We defined relative
443 dominance (D) of a climatic component (COMP) of NEE (e.g. NEE<sup>TEMP</sup>) as the mean ($\mu_s$) variance of
444 annual anomalies ($\sigma_y^2$) of this component divided by the mean variance of NEE*:

$$
445 \qquad D^{\text{COMP}} = \frac{\mu_s(\sigma_y^2(NEE_{s,y}^{\text{COMP}}))}{\mu_s(\sigma_y^2(NEE_{s,y}^{*}))} \qquad \text{EQ (6)}
$$

446 To illustrate how this relative dominance changes systematically with spatial scale we aggregated
447 NEE components successively to coarser spatial resolutions starting at 0.5° (~54.000 grid cells) and
448 ending with 'global'(1 grid cell at 360 degrees resolution) and recomputed relative dominance for
449 each spatial resolution. In total 25 levels of spatial resolution were used: 0.5, 1, 1.5, 2.5, 3, 4, 4.5, 5,
450 6, 7.5, 9, 10, 12, 15, 18, 20, 22.5, 30, 36, 45, 60, 90, 180, 360 degrees.

451 These computations were carried out for each ensemble member separately and the mean across
452 ensemble members ($\mu_e$) was plotted for each spatial resolution as dots connected with a line. The
453 uncertainty reflected by the spread of ensemble members ($\sigma_e$) was plotted as light shaded area. This
454 uncertainty is dominated by uncertainty of the mean relative dominance and not by uncertainty on
455 the systematic change with spatial aggregation. To visualize that we provided a dark shaded area in
456 the plots which represent the uncertainty on the 'shape of the curve' (U in Equation 7). This is based
457 on the standard deviation across ensemble members after subtracting the mean relative dominance
458 over all spatial resolutions (I in Equation 7) for each ensemble member (Equation 7). While Fig.2e,f
459 shows the effect of shifting relative dominance of NEE<sup>WAI</sup> vs NEE<sup>TEMP</sup> with spatial resolution
460 considering the entire global vegetated area, we repeated this analysis for different biomes (see SI-4)
461 by considering only grid cells belonging to a specific biome.

$$462 \qquad U_l = \sigma_e(D_{l,e} - \mu_l(D_{l,e})) \qquad \text{EQ (7)}$$

463 **Covariance of temperature and water availability components of GPP and TER (Fig.4).** We
464 computed the correlation coefficient and covariance between GPP and TER components (e.g. GPP<sup>TEMP</sup>
465 vs. TER<sup>TEMP</sup>) for each grid cell and ensemble member. The covariance terms were normalized to the
466 mean variance of NEE* (Equation 8). Fig. 4 shows the mean across the ensemble members (μ<sub>e</sub>) for
467 FLUXCOM, and EDF 3 the mean for the TRENDY ensemble. EDF 4 shows latitudinal patterns of the
468 spread among ensemble members (σ<sub>e</sub>) for FLUXCOM and TRENDY. The robustness of FLUXCOM
469 results with respect to different NEE flux partitioning methods is assessed in SI-6.

$$470 \qquad normalized \ COV_s(GPP_{s,y}^{COMP}, TER_{s,y}^{COMP}) = \frac{COV_y(GPP_{s,y}^{COMP}, TER_{s,y}^{COMP})}{\mu_s(\sigma_y^2(NEE_{s,y}^*))} \qquad \text{EQ (8)}$$

471 **Comparison with atmospherically based data (EDF 1).** We used three data sources of
472 atmospherically based net CO₂ flux exchange. The first is based on the annually resolved Global
473 Carbon Budget (GCP)¹³, which uses measurements of atmospheric CO₂ growth rate and estimates of
474 fossil fuel emissions, ocean uptake, and land use change emissions to derive the global land flux as a
475 residual. The second is based on the Jena CarboScope atmospheric transport inversion²⁴ (Jena
476 Inversion, version s81_3.7) covering the full time period of the study. The third is an ensemble of 10
477 atmospheric inversions¹⁹ used for the REgional Carbon Cycle Assessment and Processes (RECCAP)
478 activity covering the period 1990-2012, with each inversion covering a different time period. Four
479 versions of the Jena Inversion have been removed from the original 14 member RECCAP ensemble to
480 make it an independent assessment. We used globally integrated net land CO₂ flux estimates from
481 the three data sources to assess globally integrated NEE IAV of FLUXCOM and TRENDY. For the Jena
482 and RECCAP inversions, we additionally calculated the integrated net land CO₂ flux for areas north
483 and south of 30°N. All time series were detrended. For RECCAP inversions we calculated the median
484 estimate of the available inversion estimates per year. All time series were normalized by the
485 standard deviation of the respective globally integrated annual net land CO₂ flux.

486 **References**

487 31 Breiman, L. Random forests. Machine Learning 45, 5-32, doi:10.1023/a:1010933404324
488 (2001).
489 32 Haykin, S. Neural Networks: A Comprehensive Foundation. (Prentice Hall PTR, 1998).
490 33 Friedman, J. H. Multivariate Adaptive Regression Splines. Annals of Statistics 19, 1-67,
491 doi:10.1214/aos/1176347963 (1991).
492 34 Jung, M. & Zscheischler, J. A guided hybrid genetic algorithm for feature selection with
493 expensive cost functions. Procedia Computer Science 18, 2337-2346 (2013).
494 35 Jung, M. et al. Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and
495 sensible heat derived from eddy covariance, satellite, and meteorological observations.
496 Journal of Geophysical Research - Biogeosciences 116, G00J07, doi:10.1029/2010JG001566
497 (2011).
498 36 Jung, M. & FLUXCOM team FLUXCOM (RS+METEO) Global Land Carbon Fluxes using
499 CRUNCEP climate data. FLUXCOM Data Portal
500 http://dx.doi.org/10.17871/FLUXCOM_RS_METEO_CRUNCEPv6_1980_2013_v1 (2016)
501
502 37 Friedl, M. A. et al. MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets. *Remote Sensing of Environment* **114**, 168-182, doi:http://dx.doi.org/10.1016/j.rse.2009.08.016 (2010).
503
504
505 38 New, M., Hulme, M. & Jones, P. Representing twentieth-century space-time climate variability. Part II: Development of 1901-96 monthly grids of terrestrial surface climate. *Journal of Climate* **13**, 2217-2238, doi:10.1175/1520-0442(2000)013<2217:rpcstc>2.0.co;2 (2000).
506
507
508
509 39 Wang, Y. P., Law, R. M. & Pak, B. A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere. *Biogeosciences* **7**, 2261-2282, doi:10.5194/bg-7-2261-2010 (2010).
510
511
512 40 Jain, A. K., Meiyappan, P., Song, Y. & House, J. I. CO2 emissions from land-use change affected more by nitrogen cycle, than by the choice of land-cover data. *Global Change Biology* **19**, 2893-2906, doi:10.1111/gcb.12207 (2013).
513
514
515 41 Sitch, S. et al. Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. *Global Change Biology* **9**, 161-185 (2003).
516
517 42 Smith, B., Prentice, I. C. & Sykes, M. T. Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within European climate space. *Global Ecology and Biogeography* **10**, 621-637, doi:10.1046/j.1466-822X.2001.t01-1-00256.x (2001).
518
519
520
521 43 Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Global Biogeochemical Cycles* **19**, GB1015, doi:10.1029/2003GB002199 (2005).
522
523
524 44 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. *Journal of Land Use Science* **8**, 104-122, doi:10.1080/1747423X.2011.628705 (2013).
525
526
527
528

529 **Data availability.** The FLUXCOM data that support the findings of this study are available from the
530 Data Portal of the Max Planck Institute for Biogeochemistry (https://www.bgc-
531 jena.mpg.de/geodb/projects/Home.php) with the identifier
532 doi:10.17871/FLUXCOM_RS_METEO_CRUNCEPv6_1980_2013_v1. The TRENDY v3 data that support
533 the findings of this study are available from Stephen Sitch (S.A.Sitch@exeter.ac.uk) upon reasonable
534 request. Source data of Fig.1 a-d, Fig, 2 e-f, and Fig. 3 a-f are additionally provided as Excel
535 spreadsheets with the paper.
536

536

537 Extended Data Figure Legends

538 **Extended Data Figure 1: Global patterns of NEE IAV for FLUXCOM (left) and TRENDY (right).** *Maps* 539 of NEE IAV magnitude (mean of ensemble members, a, b) defined as standard deviation of annual 540 NEE normalized by the mean standard deviation (values above 1 indicate above average IAV). Dashed 541 lines separate areas north and south of 30°N. Time series of integrated NEE over broad latitudinal 542 bands (c-f) or global (g,h) for 1980-2013 normalized by the standard deviation of globally integrated 543 NEE. Black lines show the mean of FLUXCOM or TRENDY ensemble members and the shaded area 544 refers to the ensemble spread (1 s.d.). Independent estimates from the Global Carbon Project (GCP), 545 the Jena Inversion, and the Regional Carbon Cycle Assessment and Processes (RECCAP) inversions (see 546 Methods) are presented with coloured lines (see legend); correlation coefficients with those are given 547 in the same colour. See SI-1 for further cross-consistency analysis.
548 Extended Data Figure 2: Local vs global dominance of NEE<sup>TEMP</sup> vs NEE<sup>WAI</sup> for FLUXCOM and TRENDY ensemble members. Dots show individual ensemble members and the crosses show ensemble means with one standard deviation. Plotted is the difference of local NEE<sup>WAI</sup> and NEE<sup>TEMP</sup> dominance (difference of blue and green most left data point in Fig.2 e,f, in main article) against the difference of global NEE<sup>WAI</sup> and NEE<sup>TEMP</sup> dominance (difference of blue and green most right data point in Fig.2 e,f, in main article). The majority of ensemble members as well as ensemble means fall in the lower right quadrant meaning an overall agreement that NEE<sup>WAI</sup> dominates at individual grid cells ('locally') but NEE<sup>TEMP</sup> the globally integrated flux anomaly ('global').
549
550
551
552
553
554
555

556 Extended Data Figure 3: Spatial patterns of covariance and correlation of WAI and TEMP driven GPP and TER IAV for TRENDY models. Maps of the covariance of annual anomalies (see Equation 8 in Methods) of GPP and TER climatic components show large compensation effects (positive covariance) for WAI (a) but nearly no covariance for TEMP (c). Correlations between GPP<sup>WAI</sup> and TER<sup>WAI</sup> are large and ubiquitous positive (b) while correlations among GPP<sup>TEMP</sup> and TER<sup>TEMP</sup> are weaker with a distinct spatial pattern of negative correlations in hot regions (d). All results refer to the mean of all FLUXCOM ensemble members. See Fig.4 for equivalent FLUXCOM results, and EDF 4 for uncertainties.
557
558
559
560
561
562

563 Extended Data Figure 4: Ensemble spread of covariation between TEMP and WAI components of GPP and TER for FLUXCOM and TRENDY. Plots show mean covariance (left) and correlation (right) between GPP<sup>TEMP</sup> and TER<sup>TEMP</sup> and GPP<sup>WAI</sup> and TER<sup>WAI</sup> for latitudinal bins of 5° for individual ensemble members (thin dotted lines) and ensemble mean (thick solid line with shaded area for 1 s.d.). Despite uncertain magnitudes of GPP<sup>TEMP</sup> and TER<sup>TEMP</sup> correlation (large green shaded area in right panels) their covariance is negligible (small shaded green area in left panels). In comparison, there is large positive covariance of GPP<sup>WAI</sup> and TER<sup>WAI</sup> but its magnitude differs substantially among ensemble members (large blue shaded area in left panels).
564
565
566
567
568
569
570

571
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>