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

HAL
open science

# Lower land-use emissions responsible for increased net land carbon sink during the slow warming period

Shilong Piao, Mengtian Huang, Zhuo Liu, Xuhui Wang, Philippe Ciais, Josep

G. Canadell, Kai Wang, Ana Bastos, Pierre Friedlingstein, Richard A.

Houghton, et al.

► To cite this version:

Shilong Piao, Mengtian Huang, Zhuo Liu, Xuhui Wang, Philippe Ciais, et al. Lower land-use emissions responsible for increased net land carbon sink during the slow warming period. Nature Geoscience, 2018, 11, pp.739-743. 10.1038/s41561-018-0204-7. insu-03721862

HAL Id: insu-03721862

https://insu.hal.science/insu-03721862v1

Submitted on 17 Jun 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**ORE** Open Research Exeter

**TITLE**

Lower land-use emissions responsible for increased net land carbon sink during the slow warming period

**AUTHORS**

Piao, S; Huang, M; Liu, Z; et al.

**JOURNAL**

Nature Geoscience

**DEPOSITED IN ORE**

25 October 2018

This version available at

http://hdl.handle.net/10871/34446
1 **Lower land use emissions increased net land carbon sink during warming hiatus period**

2 Shilong Piao¹,²,³, Zhuo Liu¹, Mengtian Huang¹, Xuhui Wang¹,⁴, Philippe Ciais⁴, Josep G.

3 Canadell⁵, Kai Wang¹, Ana Bastos⁴, Pierre Friedlingstein⁶, Richard A. Houghton⁷, Corinne Le

4 Quéré⁸, Yongwen Liu¹, Ranga B. Myneni⁹, Shushi Peng¹, Julia Pongratz¹⁰, Stephen Sitch¹¹,

5 Tao Yan¹, Yilong Wang⁴, Tao Wang²,³, Zaichun Zhu¹, Donghai Wu¹

6

7 ¹ Sino-French Institute for Earth System Science, College of Urban and Environmental

8 Sciences, Peking University, Beijing 100871, China

9 ² Key Laboratory of Alpine Ecology and Biodiversity, Institute of Tibetan Plateau Research,

10 Chinese Academy of Sciences, Beijing 100085, China

11 ³ Center for Excellence in Tibetan Earth Science, Chinese Academy of Sciences, Beijing

12 100085, China

13 ⁴ Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ,

14 Gif-sur-Yvette 91191, France

15 ⁵ Carbon Project, CSIRO Oceans and Atmosphere, Canberra, Australian Capital Territory

16 2601, Australia

17 ⁶ College of Engineering, Mathematics and Physical Sciences, University of Exeter, North

18 Park Road, Exeter EX4 4QF, UK

19 ⁷ Woods Hole Research Center, 149 Woods Hole Road, Falmouth, Massachusetts 02540, USA

20 ⁸ Tyndall Centre for Climate Change Research, University of East Anglia, Norwich Research

21 Park, Norwich, NR4 7TJ, UK

22 ⁹ Department of Earth and Environment, Boston University, Boston, Massachusetts 02215,

23 USA

24 ¹⁰ Max Planck Institute for Meteorology, Bundesstraße 53, 20146 Hamburg, Germany

25 ¹¹ College of Life and Environmental Sciences, University of Exeter, North Park Road, Exeter
26 EX4 4RJ, UK

27

28
The terrestrial carbon sink has shown an acceleration after 1998, coincident with the warming hiatus¹,². However, different mechanisms were proposed¹,². Here we analyse recent change in the net land carbon sink (NLS) and its driving factors using atmospheric inversions results³,⁴ and terrestrial carbon models. We show that the linear trend of NLS during 1998-2012 (0.17±0.05 PgC yr⁻²) is three times larger than during 1980-1998 (0.05±0.05 PgC yr⁻²). This NLS intensification cannot be explained by CO₂ fertilization (0.02±0.11 PgC yr⁻²) and climate change (-0.03±0.15 PgC yr⁻²) according to terrestrial model simulation⁵,⁶. Thus, we looked more into the contribution of changes in land use emissions (E_LUC) estimated from the bookkeeping model of Houghton et al.⁷ showing decreasing E_LUC as the dominant driver (73%) of the intensification of NLS during 1998-2012. This reduction of land-use change emissions is due to both decreased tropical forest area loss and increased afforestation in northern temperate regions. Calculating E_LUC with the inversion-based estimate shows consistently reduced E_LUC, while another bookkeeping model⁸ did not reproduce such change probably due to missing the signal of reduced tropical deforestation. These results highlight the importance of better constraining emissions from land use change to understand recent trends in land carbon sinks.
Coincident with the warming hiatus of 1998-2012⁹⁻¹¹, the vegetation greening trend observed from several satellite products stalled after 1998 in most regions¹²⁻¹⁶ while the global land carbon sink has continued to increase¹,². Keenan et al.¹ and Ballantyne et al.² analysed this signal from the residual terrestrial carbon sink (RLS) calculated by difference between emissions from fossil fuel and land use, and ocean uptake and atmospheric CO₂ growth rate. The mechanisms behind the recent increase in RLS were inconsistent between the two studies. Keenan et al.¹ suggest increasing photosynthesis and decreased respiration, whereas Ballantyne et al.² suggest decreasing photosynthesis and thus reduced respiration being the only mechanism through which RLS increased during the hiatus. Furthermore, the seasonal and spatial patterns of changes in land carbon sink do not match with those of temperature changes¹⁷. Of note is the fact that systematic errors in land use emissions ⁷ directly transfer as bias of RLS⁵,¹⁸. Thus, instead of RLS, we revisit changes in the net land carbon sink (NLS) including land use emissions and its driving factors using atmospheric inversions and land carbon models.

The NLS estimated from the two inversions (see Methods) and from the global CO₂ budget¹⁹ show a three-times faster increase after 1998 (0.17±0.05 PgC yr⁻², mean ± 1 standard error) than in the decade before (0.05±0.05 PgC yr⁻²) (Fig. 1 and Supplementary Table 1, see Methods). The limit of 1998 is the one used by IPCC²⁰ and previous carbon cycle study²¹ as the beginning of the hiatus, but using 2001 or 2002 as the starting year of the warming hiatus yields similar results (Supplementary Table 2). The enlarging positive trend in NLS after 1998 (i.e. NLS intensification) is also found on a 5-years moving window (Supplementary Fig. 1) and in different inversion versions with more atmospheric CO₂ measurement sites but for shorter period (Supplementary Table 3 and Fig. 2).

NLS can be decomposed as the sum of three components, net primary productivity (NPP), heterotrophic respiration and fires in natural ecosystems (HR+F) and net carbon
emissions from land use change ($E_{LUC}$). The fraction of fire emissions that happens during land use change, known as deforestation fires, is included in $E_{LUC}$, while carbon emission from fossil fuels for land management is not included in $E_{LUC}$. To explain why NLS increased faster after 1998, we consider three mechanisms: (M1) NPP increased faster than before, forcing a sink intensification; (M2) heterotrophic respiration and fires (HR+F) increased at a slower rate than before or declined, consistent with slower warming rates; (M3) $E_{LUC}$ emissions decreased²².

**Trends in NPP** For the first mechanism, we analysed NPP changes over the past 30 years using the dynamic vegetation models (DGVM) from the TRENDY project and satellite observation-based NPP from Smith et al.¹² (hereafter SM16, see Methods). As shown in Figure 1B and 1D, both satellite-derived NPP and modelled NPP showed significant positive trends (an indication of enhanced carbon assimilation) before 1998 (SM16: 0.12±0.03 PgC yr⁻², P < 0.01; DGVMs mean: 0.15±0.04 PgC yr⁻², P < 0.01). After 1998, however, the satellite-based NPP shows a significantly (P<0.05) smaller positive trend (0.04 ± 0.04 PgC yr⁻², P > 0.05) than before. By comparison, four of the eight DGVMs do not show deceleration of NPP (i.e. reduced trend of NPP) after 1998, with trend change of NPP ranging from -0.08±0.05 PgC yr⁻² (P < 0.05) to 0.11±0.06 PgC yr⁻² (P < 0.01) (Supplementary Fig. 3). On average, the DGVMs show almost no change of NPP trend (-0.001±0.067 PgC yr⁻², P > 0.1) between the period before 1998 and that after 1998 (Fig. 2), and can thus barely explain (<1%) the intensification of NLS after 1998. A recent commentary²³ suggested that the disagreement of NPP trends between SM16 and DGVM is likely due to the underestimate of the CO₂ fertilization effect on satellite-based NPP. However, continued increase of CO₂ concentration over past three decades may not explain the intensification of NLS after 1998. The leaf area index (LAI) derived from GIMMS satellite products stalled in the recent period 1998-2012, which is not captured by DGVMs (Supplementary Fig. 4). This overestimate of
the LAI trend in the period after 1998 suggests that DGVMs may under-estimate the deceleration of NPP in the recent decade captured in SM16. Therefore, the forcing from NPP change alone cannot explain why NLS intensified.

**Trends in HR and natural fire** To analyse the second mechanism (M2) we analysed changes in HR based on the same DGVM results⁵,⁶. As shown in Fig. 2 and Supplementary Table 1, a reduction in the positive trend of HR (i.e. a deceleration of carbon emission from HR) in simulations where models were driven by changing CO₂ and climate was found by most DGVMs, with six out of the eight models showing a reduced trend of HR after 1998 ranging from -0.06±0.03 PgC yr⁻² (P < 0.01) to 0.06±0.08 PgC yr⁻² (P > 0.05). The small deceleration of HR (-0.01±0.04 PgC yr⁻², P > 0.05), however, accounts for less than 9% (-47% - 49%) of the observed intensification of NLS. According to factorial DGVM simulations, the effect of climate change alone (see Methods) did cause a significant deceleration of HR in the period 1998-2012 (-0.04±0.05 PgC yr⁻¹, P > 0.05) compared to the period 1980-1998 (Fig. 2), consistent with a slower warming rate between 1998 and 2012. However, the climate driven HR deceleration (i.e. deceleration in carbon emission) is also paralleled by a NPP deceleration (i.e. deceleration in carbon uptake) due to climate change alone in the DGVM models (-0.06±0.10 PgC yr⁻², P > 0.05; Fig. 2). This indicates that the NLS intensification during 1998-2012 cannot be attributed to climate change alone in the DGVM models. The simulation results of these models further show that rising atmospheric CO₂ can only explain 19% of the NLS intensification (Fig. 2), and that the combinations of CO₂ and climate change cancel each other. These results suggest that mechanisms other than CO₂ fertilization and climate change are responsible for the observed intensification of the NLS.

Besides HR, a reduction in natural fire emission could be another cause the intensification in the NLS. Accounting natural fires at global scale remains challenging, because satellite-based burn area cannot readily distinguish natural fires from other causes²⁴,²⁵.
Therefore, we analysed trends in fire simulated by four TRENDYv2 DGVMs, which considered wild fire processes. The models exhibited large differences in the change of fire emissions trend during the two periods (CLM4.5: -0.052±0.020 PgC yr⁻¹, P < 0.01; LPJ: 0.004±0.009 PgC yr⁻¹, P > 0.05; VISIT: 0.007±0.018 PgC yr⁻¹, P > 0.05; LPJ-GUESS: 0.013±0.024 PgC yr⁻¹, P > 0.05) (Supplementary Fig. 5). However, even considering the full model range of trend estimates, the natural fire emission probably contributes negatively to intensification of NLS (-6%±25%).

**Trends in net carbon emission from land use change** Over the last thirty years there has been a slow-down of forest losses²⁶⁻³⁰. According to the latest Forest Resources Assessment (FRA 2015) by Food and Agriculture Organization of the United Nations³¹, the annual rate of net forest loss decreased from 7.27 M ha yr⁻¹ in the 1990s to 3.99 M ha yr⁻¹ in the 2000s, primarily owing to less logging in tropical regions and increased plantations in northern temperate lands (Supplementary Table 4 and Fig. 6). Therefore, the NLS intensification can also reflect decreased E_LUC during 1998-2012.

We estimated E_LUC using the latest version of the bookkeeping model from Houghton et al.⁷ (hereafter BK), which was widely used and adopted by the Global Carbon Project in accounting annual global carbon budget³². The global E_LUC is a source of 1.13 PgC yr⁻¹, which is found mostly in tropical regions (1.31 PgC yr⁻¹), primarily Southeast Asia (0.54 PgC yr⁻¹), South America (0.38 PgC yr⁻¹) and Africa (0.38 PgC yr⁻¹) (Supplementary Fig. 7a). Tropical regions are found to be the largest contributor to global E_LUC emissions, followed by the Southern Hemisphere temperate regions as a slight source (1% of global E_LUC) (Supplementary Fig. 7a). We then compared the linear trend of E_LUC over the globe between 1980-1998 and 1998-2012. The deceleration of E_LUC contributes to a trend change of 0.09±0.01 PgC yr⁻² (P < 0.01) (Fig. 3), explaining 73% of NLS intensification. This result suggests that the faster increase of NLS after 1998 is primarily explained by decreasing E_LUC.
As shown in Fig. 3, the deceleration in global E<sub>LUC</sub> between 1980-1998 and 1998-2012 is attributed to tropical regions, where a decline of -0.08±0.01 PgC yr<sup>-2</sup> (P < 0.01) in E<sub>LUC</sub> trend is found (about 92% of the total decrease in global E<sub>LUC</sub> trend). The decline was largely in Southeast Asia (-0.05±0.01 PgC yr<sup>-2</sup>, P < 0.01) and South America (-0.016±0.004 PgC yr<sup>-2</sup>, P < 0.01) (Fig. 3), where the annual rate of net forest loss declined during the 2000s compared with 1990s<sup>31</sup>. For example, the rate of net forest loss in South America decreased from 4 M ha yr<sup>-1</sup> during the 1990s to 3.87 M ha yr<sup>-1</sup> during the 2000s, whereas the net loss rate in Southeast Asia during the 2000s (0.64 M ha yr<sup>-1</sup>) was only 30% of that during the 1990s (2.11 M ha yr<sup>-1</sup>) (Supplementary Fig. 6 and Table 4). For NH temperate regions, E<sub>LUC</sub> was found to decelerate between the two periods, with a linear trend of -0.010±0.001 PgC yr<sup>-2</sup> after 1998 (P < 0.01; about 11% of the total decrease in global E<sub>LUC</sub> trend). Temperate North America accounted for the largest fraction (89%; -0.009±0.006 PgC yr<sup>-2</sup>, P < 0.01) of decreasing E<sub>LUC</sub> in the northern temperate zone, mainly due to the fact that the forest area decrease of -0.35 M ha yr<sup>-1</sup> in the 1990s was reversed to an increase of 0.22 M ha yr<sup>-1</sup> after 2000<sup>31</sup> (Supplementary Fig. 6 and Table 4).

In addition to BK based on FAO/FRA land use areas and regional carbon response curves to land use change<sup>18</sup>, we also explored E<sub>LUC</sub> estimates with two other methods, which are the bookkeeping model of Hansis et al.<sup>8</sup> (hereafter BKH) based on Land Use Harmonization (LUH) data from 1500 to 2004<sup>33</sup> and the Global Carbon Project update from 2005 to 2012<sup>5</sup> (see Methods), and E<sub>LUC</sub> estimated by forming the difference between the net land-atmosphere CO<sub>2</sub> flux from atmospheric inversions and the fraction of this flux attributed to natural ecosystems simulated under the TRENDY S2 DGVM simulation (hereafter E<sub>Inversion-LF-DGVMs(S2)</sub>, see Methods). Globally, the change in trend of global E<sub>LUC</sub> after 1998 by E<sub>Inversion-LF-DGVMs(S2)</sub> (-0.07±0.05 PgC yr<sup>-2</sup>, P < 0.05) was similar to that by BK, but BKH estimated little change in trend of E<sub>LUC</sub> (-0.01±0.01 PgC yr<sup>-2</sup>, P > 0.05) for the same period.
The lack of trend change by BKH may come from uncertainties in land cover input dataset. Important differences between the land use input used in BK, which is directly based on FAO/FRA, and the harmonized land use dataset by Hurtt et al.³³ used in BKH are assumptions on shifting cultivation in the tropics and additional assumptions introduced in the latter dataset to make the country-level FAO/FRA data spatially explicit. Forest cover changes are not explicitly indicated by the harmonized land use dataset but deduced from changes in agricultural areas and thus can differ largely from forest inventory data both in magnitude and in trends (Supplementary Fig. 8). For example, The BKH estimated E<sub>LUC</sub> over South America exhibited positive change (0.007±0.008 PgC yr<sup>-2</sup>, P > 0.05) during the warming hiatus period, which is in contrast to forest survey data suggesting a reduced rate of deforestation in 2000s³¹.

The shift of land cover dataset in 2004 is also a potential issue making BKH more uncertain in estimating change in E<sub>LUC</sub> trend during the recent decade. The general consensus between BK and E<sub>Inversion-LF-DGVMs(S2)</sub> in estimating change of E<sub>LUC</sub> trend globally and over South America suggests the potential of utilizing this new method in estimating E<sub>LUC</sub>. However, it also differs from BK in estimating trend change of E<sub>LUC</sub> at regional scale, for example, over Africa (-0.002±0.001 PgC yr<sup>-2</sup>, P < 0.05 by BK vs. 0.04±0.03 PgC yr<sup>-2</sup>, P < 0.05 by E<sub>Inversion-LF-DGVMs(S2)</sub>; Supplementary Fig. 7b). The lack of atmospheric CO<sub>2</sub> observations over Africa can be a large source of uncertainties in atmospheric inversion, as indicated by the large error bars in regional E<sub>LUC</sub> estimates (Supplementary Fig. 7b). The uncertainties in land carbon models⁶ are also propagated in E<sub>Inversion-LF-DGVMs(S2)</sub>.

In summary, our results confirm an intensification in the NLS during the warming hiatus, (1998-2012) as compared to the preceding period (1980-1998). Using different approaches, we found that a number of drivers were responsible for the enhanced rate of the NLS. The decreasing trend in net carbon emissions from land use change was the dominant cause during warming hiatus period. The decreasing emissions from land use change were not driven by a
lower rate of warming during this period, but by reduced deforestation in the tropics and
increased afforestation in NH temperate regions. Consistent with Keenan et al.¹, we found a
lower positive trend of HR due to a lower rate of warming during the second period. But
contrary to them, our analysis, based on an ensemble of DGVMs under different scenarios
instead of a semi-empirical model¹, shows little effect of HR trends on the NLS, mainly
because of the compensating effects of CO₂ fertilization (increasing carbon emissions from
HR through higher input) and climate change (decreasing carbon emissions from HR). Note
that large uncertainties still remain with estimates of carbon flux from land use change and its
trend over the last thirty years, particularly in East Asia, South America, Africa and Europe.
Reducing this uncertainty is a top priority for future work to more accurately predict the future
evolution of the global carbon cycle and its feedback to climate change. To this end, detailed
information on LULCC transitions²⁸,³⁴ with high spatio-temporal resolution, and on carbon
response functions to these transitions³⁰,³⁵ is needed. In addition, various forms of land use
management (e.g. wood harvest, shifting cultivation, cropland management, fire management,
peatland drainage) are often inconsistently and incompletely represented in DGVMs⁵,¹⁸. A
better characterization of these critical processes is required in future studies.

213
214 Methods

215
**Satellite-based NDVI and NPP data.** The Normalized Difference Vegetation Index (NDVI),
216
which has been widely used to monitor vegetation activity, was obtained from Global
217
Inventory Modelling and Mapping Studies (GIMMS) third-generation product (NDVI<sub>3g</sub>) at a
218
resolution of 8 km×8 km from 1982 to 2015<sup>36</sup>.

219
The satellite-derived net primary productivity (NPP) was from MODIS<sup>13</sup> and a recent
220
study by Smith et al.<sup>12</sup> (SM16). For the latter, NPP was calculated based on MODIS NPP
221
algorithm<sup>13</sup>, but driven by 30-year (1982-2011) GIMMS fraction of photosynthetically active
222
radiation (FPAR) and leaf area index (LAI) data<sup>12</sup>. Further details about satellite-derived NPP
223
data can be found in Smith et al.<sup>12</sup> and Zhao & Running<sup>13</sup>. Note that the MODIS results only
224
cover the period from 2001 onwards. Therefore, we only included the MODIS results in
225
Supplementary Fig. 9 to show that the stall of NPP during warming hiatus period is not an
226
artifact from the only one long-term satellite-derived net primary productivity (NPP) data
227
from Smith et al.<sup>12</sup>.

228
**Dynamic global vegetation models (DGVMs).** An ensemble of eight dynamic global
229
vegetation models (Supplementary Table 5 from the project "Trends and drivers of the
230
regional scale sources and sinks of carbon dioxide" (TRENDY) were used to simulate the
231
carbon balance of terrestrial ecosystems during the period 1980-2012. These models provided
232
outputs of Net Biome Productivity (NBP), Net Primary Productivity (NPP) and Heterotrophic
233
Respiration (HR). Here we used NBP to reflect the magnitude of net land carbon sink (NLS,
234
NLS = NBP = NPP - HR - D, D refers to other losses of carbon due to disturbance, including
235
carbon emissions from land use change). Note that we adopted the convention that a sink of
236
CO₂ is defined as positive (removing CO₂ from the atmosphere).

237
The DGVMs were coordinated to perform three simulations (S1, S2 and S3) following
238
the TRENDY protocol<sup>6</sup>. In simulation S1, only atmospheric CO<sub>2</sub> concentration was varied. In
simulation S2, atmospheric CO₂ and climate were varied. In simulation S3, atmospheric CO₂, climate and land use were varied. The effects of rising atmospheric CO₂, climate change and land use change on NLS can then be obtained from S1, the difference between S2 and S1, and the difference between S3 and S2, respectively. All models used the same forcing datasets, of which global atmospheric CO₂ concentration was from the combination of ice core records and atmospheric observations³⁷; historical climate fields were from CRU-NCEP dataset (http://dods-extra.cea.fr/data/p529viov/cruncep/); land use data were from the Land Use Harmonization dataset³¹ based on the History Database of the Global Environment (HYDE)³⁸. All the model outputs were resampled to a spatial resolution of 0.5°×0.5° based on the nearest neighbour method.

Note that there is a large difference between TRENDYv2 and TRENDYv4 in the estimate of NLS trend before and after 1998 under S3 simulation (Supplementary Fig. 10). On average, NLS in TRENDYv2 shows a non-significant trend before 1998 and a significant increasing trend after 1998 (Supplementary Fig. 10h), which is consistent with the results from the global carbon budget and atmospheric inversions. However, in TRENDYv4, an opposite case was found (Supplementary Fig. 10h). This difference between TRENDYv2 and TRENDYv4 in simulating the observed NLS trend mainly results from the simulation of land use change rather than S2 simulation (Supplementary Fig. 10h). This not only indicates large uncertainties in the simulation of land use change (Supplementary Fig. 7), but suggests the potential effect of land use change on NLS trend. Although TRENDYv4 used an updated and improved input of land use change maps (HYDE3.2)³⁹ compared with TRENDYv2 (HYDE3.1), we did not adopt it to estimate carbon emissions from land use change given that it did not capture the trend of NLS before and after 1998. Overall, we only used TRENDY results derived from S1 and S2 simulation in our main text, and proposed a new way to estimate land use change emission by combining the results from atmospheric inversions and TRENDY models under
264 S2 simulation (see below).

265 **Global carbon budget.** To gain a better understanding of the net land carbon sink, we also
266 used data from global carbon budget coordinated by the Global Carbon Project (GCP)¹⁹. Here
267 the net land sink was inferred as a residual of fossil fuel emissions, atmospheric CO₂
268 accumulation and ocean sink, which is independent from atmospheric inversions.

269 **Atmospheric CO₂ inversion data.** Atmospheric CO₂ inversions offer a method in which CO₂
270 observation networks, transport models and a prior knowledge of fluxes are utilized to
271 estimate net land-atmosphere carbon exchange⁴⁰. This top-down approach allows us to
272 compare the magnitude of net land carbon sink (NLS) with that from bottom-up method based
273 on DGVMs. Given our long-term study period from 1980 to 2012, here we used two inversion
274 products: MACC_v15 from Chevallier et al.³ (hereafter MACC, available time period:
275 1979-2015) and JENA_S81_v3.8 from Rödenbeck et al.⁴ (hereafter JENA, available time
276 period: 1981-2014). The original spatial resolution of MACC and JENA is
277 1.875°latitude×3.75°longitude and 3.75°latitude×5°longitude, respectively.

278 It should be noted that there are differences between these two inversions in number of
279 observation sites as constraint, transport models and prior flux information⁴⁰. As
280 recommended in previous studies⁴⁰,⁴¹, a standard fossil fuel and cement production flux (FFC)
281 should be subtracted from the total posterior fluxes when comparing net land flux from
282 different CO₂ inversions. This is due to the fact that differences in prior FFC will manifest as
283 differences in the estimated natural flux³⁸. Thus, here we took the fossil fuel flux which is
284 used in GCP carbon budget as a standard and subtracted it from the total posterior fluxes for
285 both CO₂ inversions to obtain the “fossil corrected” NLS, although the global fossil fuel
286 emissions are quite consistent between the two inversions and with the GCP data
287 (Supplementary Fig. 11). Note that the FFC data used in GCP carbon budget was from the
288 Carbon Dioxide Information Analysis Center (CDIAC,
289

http://cdiac.ornl.gov/trends/emis/meth_reg.html) and energy statistics published by BP
(http://www.bp.com/en/global/corporate/about-bp.html).

290

291

**Net carbon flux from land use change (E<sub>LUC</sub>)**. We used the estimates by Houghton et al.⁷ (hereafter BK) for carbon fluxes due to land use change. In this method, ground-based measurements of carbon density are combined with land cover change data from the Forest Resource Assessment (FRA) of the Food and Agriculture Organization (FAO) using a semi-empirical bookkeeping model, in which standard growth and decomposition curves are used to track changes in carbon pools¹⁸. Using the estimate by Houghton et al.⁷ is consistent with the global carbon budget estimates provided by the Global Carbon Project⁴², but may conceal large uncertainties associated with land use change itself as well as LUC-related carbon fluxes. We therefore include in the supplemental analyses two additional approaches: The second approach is also a bookkeeping method but from Hansis et al.⁸ (hereafter BKH). Although BKH largely follows the bookkeeping method developed by Houghton et al.⁴³,⁴⁴, there are key differences between BKH and BK: BKH is spatially explicit at a resolution of 0.5°×0.5°⁸, whereas BK is constructed based on aggregated, non-spatial national and international statistics¹⁸; BKH used Land Use Harmonization dataset from 1500 to 2004³¹ and the Global Carbon Project update from 2005 to 2012 as input⁸ while BK used FAO/FRA land use change data¹⁸; other differences between BKH and BK are the accounting of successive LULCC events including their interactions in BKH and different assumptions on the allocation of agricultural land on natural vegetation⁸. Note that the data available now from Houghton et al.⁴³,⁴⁴ and Hansis et al.⁸ does not enable us to obtain the quantifiable uncertainties for trends.

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

Apart from above two bookkeeping approaches, here we developed a new way to
indirectly estimate E<sub>LUC</sub> using the difference of land carbon flux from atmospheric inversions,
the flux from lateral transport (LF) and that from DGVMs under S2 simulation (driven by

312

313
rising CO₂ and climate change, not taking into account LF) (hereafter referred to EInversion-LF-DGVMs(S2)). This approach was based on the assumption that the effect of changing atmospheric CO₂ concentration and climate are well modelled by DGVMs so that the difference between inversion fluxes (including all CO₂ sources and components), lateral carbon flux and DGVM modelled fluxes under S2 simulation equals the net source from land use and land management.

The processes of lateral carbon transport generally involve (1) the trade of food and wood products; (2) carbon export from land to ocean by rivers. In terms of the lateral carbon flux associated with food and wood trade (Supplementary Fig. 12), we first derive the annual import and export data of food and wood products from FAO statistical databases (http://www.fao.org/faostat/en/#data). Then the food and wood data are converted into dry biomass and into carbon using specific conversion factors. For food products, we adopted crop-specific coefficients (including dry matter content of harvested biomass and carbon content of harvested dry matter, see Supplementary Table 6) following Wolf et al.⁴⁵ and Kyle et al.⁴⁶. For wood products, we adopted an average wood density of 0.5 and 0.45 carbon concentration in dry biomass following Ciais et al.⁴⁷. In terms of the carbon exported from ecosystems by rivers, we included dissolved organic carbon (DOC), particulate organic carbon (POC) and dissolved inorganic carbon (DIC) from 45 major zones (MARCATS: MARgins and CATchments Segmentation) and 149 sub-units (COSCATs: Coastal Segmentation and related CATchments)⁴⁸,⁴⁹ (http://www.biogeomod.net/geomaps/, see Supplementary Table 7). Then we aggregated the riverine carbon transport into continental scale (Supplementary Fig. 13). However, it should be noted that the carbon transport data is only a rough estimate and lack temporal evolution. Besides, it is unclear whether the exported carbon by rivers is from old deposits or from current photosynthesis. In addition, time series of the carbon exports from rivers are not available. Therefore, we did not count this part in the
calculation of LF.

Note that we obtained eighteen estimates from $F_{\text{Inversion-LF-DGVMs(S2)}}$ approach, as eight DGVMs and two atmospheric inversions were considered in the analysis. All datasets from atmospheric inversions and DGVMs were first regridded into a common $0.5°×0.5°$ grid using nearest neighbor interpolation method. We also performed the same analyses by regridding all the datasets into a common $1°×1°$ or $2°×2°$ grid, and found similar results (Supplementary Fig. 14). In addition, given that BK was based on national data and not spatially explicit, we obtained latitudinal results (the bottom left in Fig. 3) by roughly aggregating northern North America, Europe and Asian Russia into boreal region, southern North America, West/Central/South Asia and East Asia to Northern Hemisphere (NH) temperate region, South America, Africa and Southeast Asia to tropics, and Oceania to Southern Hemisphere (SH) temperate region.

There is a S3 simulation of TRENDY where DGVMs are driven by the land cover dataset (LUH) in addition to change in climate and atmospheric CO$_2$. Thus, the difference of S3 and S2 simulations may also represent the model simulated emission of land use change. However, comparing the difference between S3 and S2 and $E_{\text{LUC}}$ estimated by the bookkeeping or inversion-based approach are difficult, because DGVMs do not simulate the full range of processes related to $E_{\text{LUC}}$ (not all DGVMs account for example for wood and crop harvest or shifting cultivation$^{42}$). Further, land use change emissions derived as difference between S3 and S2 differ in the terms that are included as compared to other approaches$^{48}$. Most notably, the loss of additional sink capacity is attributed to $E_{\text{LUC}}$ using S3 minus S2, while it is excluded from $E_{\text{LUC}}$ derived from bookkeeping models or the inversion-based approach. Lastly, the input land cover dataset has discontinuity issue in the recent decade and different models also have different assumption converting LUH dataset into model-specific land cover inputs, making it less reliable in estimating trend in the recent
decade. Therefore, we do not include the difference of S3 and S2 simulation by DGVMs in
this study.

**Statistical analysis.** We calculated the trend of NLS, NPP, HR, NDVI and E<sub>LUC</sub> during three study periods (1980-2012, 1980-1998, and 1998-2012) based on Linear Least Square Regression analysis, in which above five indicators were regarded as dependent variables and year as independent variable. The slope of the regression was then defined as the trend. The standard error of linear regression coefficient (slope) was defined as the uncertainty of the linear trend. Note that for the average trend of different data sources, the uncertainty of its trend was estimated as the root-mean-square of the standard error of for each data sources under the assumption that data from different datasets is independent from each other. Based on this, we obtained the change of above five indicators' trend between the second period (1998-2012) and the first period (1980-1998). The dividing year 1998 is selected according to IPCC description of the warming hiatus period<sup>20</sup>. However, the intensification of NLS and dominant contribution of E<sub>LUC</sub> will not change, if trend analyses starts from 2001/2002 after the El Nino/La Nina events at the end of 20<sup>th</sup> century (Supplementary Table 2). Note that here changes in the intensity of each component of NLS were indicated by changes in the magnitude (absolute value) of each term. In this case, a positive trend in NPP / HR, F and E<sub>LUC</sub> refers to an increase of carbon assimilation / carbon emission, and vice versa, a negative trend in NPP / HR, F and E<sub>LUC</sub> indicates a decline in carbon assimilation / carbon emission. The statistics of the change in trend for each flux was estimated using bootstrap analyses<sup>51</sup>. We first obtained probability distribution of NLS trend before and after 1998 in 500-time bootstrapping. Then the probability distribution in the change in trend for each flux was calculated based on the differences of trends among the sampling of the two probability distributions. For clarification, NLS intensification indicates increase in the trend of NLS after 1998. Similarly, acceleration/deceleration of a flux (NPP, HR, fire and E<sub>LUC</sub>) indicates
389 larger/smaller trend of the flux during 1998-2012 than that during 1980s-1998.

390 **Data availability.** The GIMMS NDVI<sub>3g</sub> datasets are available 391 at http://ecocast.arc.nasa.gov/data/pub/gimms/3g.v0/. The satellite-derived NPP dataset is 392 available on request from W. K. Smith<sup>12</sup>. The MODIS NPP dataset is available on request 393 from M. Zhao<sup>13</sup>. Net carbon flux from land use change (E<sub>LUC</sub>) estimated using the 394 bookkeeping approach is available on request from R. A. Houghton<sup>7</sup> and E. Hansis<sup>8</sup>, 395 respectively. Model outputs generated by Dynamic Global Vegetation Model (DGVM) 396 groups are available from Stephen Stich (s.a.sitch@exeter.ac.uk) or Pierre Friedlingstein 397 (p.friedlingstein@exeter.ac.uk) upon request.

398

399
400 **Reference**

1. Keenan, T. F. et al. Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake. *Nat Commun* **7**, 13428 (2016).

2. Ballantyne, A. et al. Accelerating net terrestrial carbon uptake during the warming hiatus due to reduced respiration. *Nat. Clim. Change* **7**, 148-152 (2017).

3. Chevallier, F. et al. CO2 surface fluxes at grid point scale estimated from a global 21 year reanalysis of atmospheric measurements. *J. Geophys. Res.* **115**, D21307, doi:10.1029/2010JD013887 (2010).

4. Rödenbeck, C. Estimating CO₂ sources and sinks from atmospheric mixing ratio measurements using a global inversion of atmospheric transport, Technical Report 6, Max Planck Institute for Biogeochemistry, Jena, available at: http://www.bgc-jena.mpg.de/uploads/Publications/TechnicalReports/tech_report6.pdf, 2005.

5. Le Quéré, C. et al. Global carbon budget 2013. *Earth Syst. Sci. Data* **6**, 235-263 (2014).

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

7. Houghton, R. A. & Nassikas, A. A. Global and regional fluxes of carbon from land use and land cover change 1850-2015. *Glob. Biogeochem. Cycles* **31**, doi:10.1002/2016GB005546 (2017).

8. Hansis, E., Davis, S. J. & Pongratz, J. Relevance of methodological choices for accounting of land use change carbon fluxes. *Glob. Biogeochem. Cycles* **29**, 1230-1246 (2015).

9. Easterling, D. R. & Wehner, M. F. Is the climate warming or cooling? *Geophys. Res. Lett.* **36**, L08706, doi:10.1029/2009GL037810 (2009).
423

10. Kaufmann, R. K., Kauppi, H., Mann, M. L. & Stock, J. H. Reconciling anthropogenic climate change with observed temperature 1998-2008. *Proc. Natl Acad. Sci. USA.* **108**, 11790-11793 (2011).

424

425

426

11. Cohen, J. L., Furtado, J. C., Barlow, M. & Alexeev, V. A. Asymmetric seasonal temperature trends. *Geophys. Res. Lett.* **39**, L04705, doi:10.1029/2011GL050582 (2012).

427

428

12. Smith, W. K. et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO₂ fertilization. *Nat. Clim. Change* **6**, 306-310 (2016).

429

430

13. Zhao, M. & Running, S. W. Drought-induced reduction in global terrestrial net primary production from 2000 through 2009. *Science* **329**, 940-943 (2010).

431

432

14. Jong, R., Verbesselt, J., Schaepman, M. E. & Bruin, S. Trend changes in global greening and browning: contribution of short-term trends to longer-term change. *Glob. Change Biol.* **18**, 642-655 (2012).

433

434

435

15. Mohammat, A. et al. Drought and spring cooling induced recent decrease in vegetation growth in Inner Asia. *Agric. For. Meteorol.* **178**: 21-30 (2013).

436

437

16. Kong, D., Zhang, Q., Singh, V. P. & Shi, P. Seasonal vegetation response to climate change in the Northern Hemisphere (1982–2013). *Glob. Planet. Change* **148**, 1-8 (2017).

438

439

17. Zhu, Z. et al. The accelerating land carbon uptake of the 2000s may not be driven predominantly by the warming hiatus. *Geophys Res Lett.* **45**: 1402-1409 (2018).

440

441

18. Houghton, R. A. et al. Carbon emissions from land use and land-cover change. *Biogeosciences* **9**, 5125-5142 (2012).

442

443

19. Le Quéré, C. et al. Global carbon budget 2015. *Earth Syst. Sci. Data* **7**, 349-396 (2015).

444

20. Hartmann, D. L. et al. Observations: Atmosphere and Surface. in Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) 192–194 (Cambridge Univ. Press, Cambridge, 2013).

445

446

447

21. Ballantyne, A. et al. Accelerating net terrestrial carbon uptake during the warming hiatus
448

due to reduced respiration. *Nature Climate Change* **7**, 148-152 (2017).

449

22. Grassi, G. et al. The key role of forests in meeting climate targets requires science for credible mitigation. *Nat. Clim. Change* **7**, 220-226 (2017).

450

451

23. De Kauwe, M. G. et al. Satellite based estimates underestimate the effect of CO₂ fertilization on net primary productivity. *Nat. Clim. Change* **6**, 892-893 (2016).

452

453

24. Giglio, L., Randerson, J. T. & van der Werf, G. R. Analysis of daily, monthly, and annual burned area using the fourth-generation global fire emissions database (GFED4). *J. Geophys. Res. Biogeosci.* **118**, 317-328 (2013).

454

455

456

25. Andela, N. et al. A human-driven decline in global burned area. *Science* **356**, 1356-1362 (2017).

457

458

26. Rudel, T. K. et al. Forest transitions: towards a global understanding of land use change. *Glob Environ Change* **15**, 23-31 (2005).

459

460

27. Sánchezcuervo, A. M., Aide, T. M., Clark, M. L. & Etter, A. Land Cover Change in Colombia: Surprising Forest Recovery Trends between 2001 and 2010. *PLoS One* **7**, e43943 (2012).

461

462

463

28. Magliocca, N. R. et al. Synthesis in land change science: methodological patterns, challenges, and guidelines. *Reg Environ Change* **15**, 211-226 (2015).

464

465

29. Chazdon, R. L. et al. Carbon sequestration potential of second-growth forest regeneration in the Latin American tropics. *Sci Adv* **2**, doi: 10.1126/sciadv.1501639 (2016).

466

467

30. Poorter, L. et al. Biomass resilience of Neotropical secondary forests. *Nature* **530**, 211 (2016).

468

469

31. Food and Agriculture Organization of the United States (FAO). Global Forest Resources Assessment 2015: How are the world's forests changing? http://www.fao.org/forest-resources-assessment/en/ (2015).

470

471

472

32. Le Quéré, C. et al. Global carbon budget 2014. *Earth Syst Sci Data* **7**, 47-85 (2015).
473

33. Hurtt, G. C. et al. Harmonization of land-use scenarios for the period 1500-2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands. *Climatic change* **109**, 117-161 (2011).

474

475

476

34. Erb, K. H. et al. Bias in the attribution of forest carbon sinks. *Nat. Clim. Change* **3**, 854-856 (2013).

477

478

35. Poeplau, C. et al. Temporal dynamics of soil organic carbon after land-use change in the temperate zone-carbon response functions as a model approach. *Glob. Change Biol.* **17**, 2415-2427 (2011).

479

480

481

36. Tucker, C. J. et al. An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. *Int. J. Remote Sens.* **26**, 4485–4498 (2005).

482

483

37. Keeling, C. D. & Whorf, T. P. Atmospheric carbon dioxide record from Mauna Loa. *Trends: A Compendium of Data on Global Change. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, Oak Ridge, TN (2005).*

484

485

486

38. Klein Goldewijk, K., Beusen, A., Van Drecht, G. & De Vos, Martine. The HYDE 3.1 spatially explicit database of human-induced global land-use change over the past 12,000 years. *Global Ecol. Biogeogr.* **20**, 73-86 (2011).

487

488

489

39. Klein Goldewijk, K. A historical land use data set for the Holocene; HYDE 3.2. *EGU General Assembly Conference Abstracts* **18**, 1574 (2016).

490

491

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

492

493

41. Thompson, R. L. et al. Top-down assessment of the Asian carbon budget since the mid 1990s. *Nat. Commun.* **7**, DOI: 10.1038/ncomms10724 (2016).

494

495

42. Le Quéré, C. et al. Global carbon budget 2016. *Earth Syst Sci Data* **8**, 605 (2016).

496

43. Houghton, R. A. et al. Changes in the Carbon Content of Terrestrial Biota and Soils between 1860 and 1980: A Net Release of CO₂ to the Atmosphere. *Ecological monographs*

497
498 **53**, 235-262 (1983).

499 44. Houghton, R. A. Revised estimates of the annual net flux of carbon to the atmosphere from changes in land use and land management 1850–2000. *Tellus B* **55**, 378-390 (2003).

501 45. Wolf, J. et al. Biogenic carbon fluxes from global agricultural production and consumption. *Glob. Biogeochem. Cycles* **29**, 1617-1639 (2015).

503 46. Kyle, P. et al. GCAM 3.0 agriculture and land use: data sources and methods. Pacific Northwest National Laboratory. PNNL-21025 (2011).

505 47. Ciais, P. et al. The impact of lateral carbon fluxes on the European carbon balance. *Biogeosciences Discuss* **3**, 1529-1559 (2006).

507 48. Laruelle, G. G. et al. Global multi-scale segmentation of continental and coastal waters from the watersheds to the continental margins. *Hydrol. Earth Syst. Sc.* **17**, 2029 (2013).

509 49. Regnier, P. et al. Anthropogenic perturbation of the carbon fluxes from land to ocean. *Nature Geosci.* **6**, 597-607 (2013).

511 50. Pongratz, J., Reick, C. H., Houghton, R. & House, J. Terminology as a key uncertainty in net land use and land cover change carbon flux estimates. *Earth System Dynamics* **5**, 177 (2014).

514 51. Manly B. F. J. Randomization, bootstrap and Monte Carlo methods in biology (CRC Press, 2006).

516

517 Acknowledgements

518 This study was supported by the National Natural Science Foundation of China (41530528), the 111 Project (B14001), and the National Youth Top-notch Talent Support Program in China.

520 We thank the TRENDY modelling group for providing the model simulation data.

521

522 Author Contributions
523

S.Piao designed the study. Z.L performed the analysis. S.Piao and Z.L drafted the paper. All authors contributed to the interpretation of the results and to the text.

524

525

526

**Author Information**

527

The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to S.Piao (slpiao@pku.edu.cn).

528

529
530 **Figure legends**

531 **Figure 1 Anomalies and liner trends of global annual net land carbon sink (NLS) (a, c)**
532 and **net primary productivity (NPP) (b, d)**. Our whole study period is from 1980 to 2012,
533 and we calculated the trends of above variables for three time periods: 1980-2012, 1980-1998
534 and 1998-2012. In the left panels, positive value refers to a net carbon sink, while negative
535 value refers to a net carbon source. The shaded area in the left panels indicates data
536 uncertainty (±1σ). In the right panel, we denote significant trends (P < 0.05) with two asterisk
537 based on t test. The error bars in the right panels indicate the standard error of linear trend for
538 each dataset. In panel (d), the range of the data (minimum-maximum range) across different
539 models is given as colored vertical bars with the solid line showing the average value. Note
540 that different colors correspond to different sources of data (see Methods), which are noted in
541 the legends of each panel.

542 **Figure 2 Change in the trend of net land carbon sink (NLS), net primary productivity (NPP) and heterotrophic respiration (HR) estimated by eight Dynamic Global Vegetation Models (DGVMs) under different scenarios between 1998-2012 and 1980-1998.** For each model, the change in the trend of NLS / NPP / HR were obtained as the trend of each variable during 1998-2012 minus that during 1980-1998. Results for the effect of rising atmospheric CO₂ concentration ('CO₂'), climate change ('CLI'), and above two factors combined ('CO₂+CLI') are shown. On each box, the central line marks the median, the edges of the box correspond to the 25th and 75th percentiles, and the whiskers extend to the range of the data. The solid dot shows the average value of the model results.

551 **Figure 3 Linear trend of net carbon emission from land use change (ELUC) and change in ELUC trend between 1998-2012 and 1980-1998.** The bottom left show results at latitudinal scale, including boreal (50°N-90°N), northern temperate (23°N-50°N), tropical (23°N-23°S) and southern temperate region (23°S-60°S). The ELUC trend during each of the two periods as
555

556

557

558

559

560

561

562

well as change in E_LUC trend between two periods are obtained based on annual E_LUC from the bookkeeping method (BK, see Methods). A positive trend refers to increased E_LUC during corresponding period, while a negative trend refers to decreased E_LUC during corresponding period. The error bars indicate the uncertainty for E_LUC trend / the change in E_LUC trend. The uncertainty of the linear trend was estimated as the standard error of linear regression coefficient (slope), while the uncertainty of the change in E_LUC trend was estimated using bootstrap analyses (see Methods).
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>