# Compatible fossil fuel CO₂ emissions in the CMIP6 earth system models’ historical and shared socioeconomic pathway experiments of the twenty-first century

Spencer K. Liddicoat, Andy J. Wiltshire, Chris D. Jones, Vivek K. Arora, Victor Brovkin, Patricia Cadule, Tomohiro Hajima, David M. Lawrence, Julia Pongratz, Jörg Schwinger, et al.

► To cite this version:

Spencer K. Liddicoat, Andy J. Wiltshire, Chris D. Jones, Vivek K. Arora, Victor Brovkin, et al. Compatible fossil fuel CO₂ emissions in the CMIP6 earth system models’ historical and shared socioeconomic pathway experiments of the twenty-first century. Journal of Climate, 2021, 34, pp.2853-2875. 10.1175/JCLI-D-19-0991.1 . insu-03668381

HAL Id: insu-03668381

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

Submitted on 15 May 2022

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Compatible Fossil Fuel CO₂ Emissions in the CMIP6 Earth System Models’ Historical and Shared Socioeconomic Pathway Experiments of the Twenty-First Century

SPENCER K. LIDDICOAT,ª ANDY J. WILTSHIRE,ª CHRIS D. JONES,ª VIVEK K. ARORA,ᵇ VICTOR BROVKIN,ᶜ,ᵈ

PATRICIA CADULE,ᵉ TOMOHIRO HAJIMA,ᵉ DAVID M. LAWRENCE,ᵉ JULIA PONGRATZ,ᶜ, JÖRG SCHWINGER,<fnref n="1" />

ROLAND SÉFÉRIAN,ⱼ JERRY F. TJIPUTRA,<fnref n="1" /> AND TILO ZIEHN<fnref n="1" />

(Manuscript received 25 December 2019, in final form 4 December 2020)

**ABSTRACT:** We present the compatible CO₂ emissions from fossil fuel (FF) burning and industry, calculated from the historical and Shared Socioeconomic Pathway (SSP) experiments of nine Earth system models (ESMs) participating in phase 6 of the Coupled Model Intercomparison Project (CMIP6). The multimodel mean FF emissions match the historical record well and are close to the data-based estimate of cumulative emissions (394 ± 59 GtC vs 400 ± 20 GtC, respectively). Only two models fall inside the observed uncertainty range; while two exceed the upper bound, five fall slightly below the lower bound, due primarily to the plateau in CO₂ concentration in the 1940s. The ESMs' diagnosed FF emission rates are consistent with those generated by the integrated assessment models (IAMs) from which the SSPs' CO₂ concentration pathways were constructed; the simpler IAMs' emissions lie within the ESMs' spread for seven of the eight SSP experiments, the other being only marginally lower, providing confidence in the relationship between the IAMs' FF emission rates and concentration pathways. The ESMs require fossil fuel emissions to reduce to zero and subsequently become negative in SSP1-1.9, SSP1-2.6, SSP4-3.4, and SSP5-3.4over. We also present the ocean and land carbon cycle responses of the ESMs in the historical and SSP scenarios. The models' ocean carbon cycle responses are in close agreement, but there is considerable spread in their land carbon cycle responses. Land-use and land-cover change emissions have a strong influence over the magnitude of diagnosed fossil fuel emissions, with the suggestion of an inverse relationship between the two.

**KEYWORDS:** Carbon cycle; Carbon dioxide; Climate models; General circulation models; Model comparison

## 1. Introduction

The release of carbon dioxide and other greenhouse gases since the onset of the industrial era has warmed the planet by around<fnref n="1" /> <fnref n="1" />°C (Stocker et al. 2013). The changing climate is already affecting natural terrestrial and marine ecosystems, thawing permafrost, and altering the hydrological cycle, with implications for drought, flooding, and wildfires. Human systems are being impacted increasingly, with reductions in crop yields, disruption

to infrastructure, water supply and food production, and increased mortality due to high and low temperature extremes (IPCC 2014). With global emissions of CO₂ continuing to rise (Friedlingstein et al. 2019), the 2015 Paris Agreement of the United Nations Framework Convention on Climate Change (UNFCCC) saw 196 governments around the world agree to a commitment to constrain greenhouse gas (GHG) emissions to a level that would limit long-term warming, relative to the preindustrial level, to well below 2°C, ideally to below 1.5°C, with emissions peaking as soon as possible and then reducing to zero within decades (UNFCCC 2015).

In its Fifth Assessment Report, the Intergovernmental Panel on Climate Change (Collins et al. 2013) reported for the first time that the change in global mean temperature relative to the preindustrial is proportional to the cumulative total CO₂ emissions released, albeit with significant uncertainty. Whether warming is limited to 1.5° or 2°C above preindustrial levels is therefore greatly dependent on how CO₂ emissions proceed over the coming century. How society evolves in the future, how the global population will change, and the extent to

<footnote>ª Met Office Hadley Centre, Exeter, United Kingdom</footnote>

<footnote>ᵇ Canadian Centre for Climate Modelling and Analysis, Environment Canada, University of Victoria, Victoria, British Columbia, Canada</footnote>

<footnote>ᶜ Max Planck Institute for Meteorology, Hamburg, Germany</footnote>

<footnote>ᵈ CEN, Universität Hamburg, Hamburg, Germany</footnote>

<footnote>ᵉ IPSL, CNRS, Sorbonne Université, Paris, France</footnote>

<footnote>f Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan</footnote>

<footnote>g Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, Colorado</footnote>

<footnote>h Ludwig-Maximilian University, Munich, Germany</footnote>

<footnote>i NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Bergen, Norway</footnote>

<footnote>j CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France</footnote>

<footnote>k CSIRO Oceans and Atmosphere, Aspendale, Victoria, Australia</footnote>

<footnote n="1">¹ 1.01°C relative to 1850–79 (Haustein et al. 2017).</footnote>

<footnote>Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0991.s1.</footnote>

<footnote>*Corresponding author:* Spencer K. Liddicoat, spencer.liddicoat@metoffice.gov.uk</footnote>
<table><caption>TABLE 1. Summary of ACCESS-ESM1.5, CanESM5, and CESM2. See supplementary information for detailed model descriptions and explanation of acronyms used in the table.</caption><thead><tr><th></th><th>ACCESS-ESM1.5</th><th>CanESM5</th><th>CESM2</th></tr></thead><tbody><tr><td>Atmosphere</td><td>UM7.3</td><td>CanAM5</td><td>CAM6</td></tr><tr><td>Resolution, vertical levels</td><td>1.875° × 1.25°, 38</td><td>T63, 2.8°, 49</td><td>0.9° × 1.25°, 32</td></tr><tr><td>Land</td><td>CABLE2.4 (CASA-CNP)</td><td>CLASS and CTEM</td><td>CLM5</td></tr><tr><td>Dynamic vegetation</td><td>No</td><td>No</td><td>No</td></tr><tr><td>Natural vegetation PFTs</td><td>4 trees, 2 grasses, 1 shrub, wetlands, tundra</td><td>5 trees, 2 grasses</td><td>8 trees, 3 shrubs, 4 grasses</td></tr><tr><td>Soil carbon pools</td><td>3</td><td>1</td><td>3</td></tr><tr><td>Litter carbon pools</td><td>3</td><td>1</td><td>4</td></tr><tr><td>Net or gross transitions</td><td>Net</td><td>Net</td><td>Net</td></tr><tr><td>Land-use PFTs</td><td>1 crop</td><td>2 crops</td><td>6 crops</td></tr><tr><td>Harvest</td><td>No</td><td>No</td><td>Crop and wood harvest</td></tr><tr><td>Grazing</td><td>No</td><td>No</td><td>No</td></tr><tr><td>Fire-carbon interactions</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>Nitrogen cycle</td><td>Yes, and phosphorus</td><td>No</td><td>Yes</td></tr><tr><td>Nitrogen limitation of carbon uptake</td><td>Yes</td><td>Downregulation of photosynthesis as CO<sub>2</sub> rises</td><td>FUN model simulates carbon cost of N uptake</td></tr><tr><td>Ocean</td><td>MOM5</td><td>NEMO 3.4.1</td><td>POP2</td></tr><tr><td>Resolution</td><td>1°</td><td>1°</td><td>1°</td></tr><tr><td>Ocean biogeochemistry</td><td>WOMBAT</td><td>CMOC</td><td>MARBL</td></tr></tbody></table>

which countries and communities will cooperate to reduce emissions from the energy sector, transport, and agriculture cannot be known with any certainty. However, integrated assessment models (IAMs), computer models of intermediate complexity, can be used to estimate how societal changes such as population growth, economic development, investment priorities, and degree of international cooperation will impact future greenhouse gas emission rates and concentrations (Gidden et al. 2019).

Phase 6 of the Coupled Model Intercomparison Project (CMIP6) is divided into several individual model intercomparison projects (MIPs), each with a particular science focus (Eyring et al. 2016), undertaking experiments with full-complexity Earth system models (ESMs) and general climate models. ScenarioMIP (O'Neill et al. 2016) is designed to compare ESMs' responses to eight new forcing scenarios out to 2100 built around five Shared Socioeconomic Pathway (SSP) scenarios of potential future social and economic development (Riahi et al. 2017).

The current generation of Earth system models typically contain fully process-based representations of the terrestrial and ocean carbon cycles, allowing uptake of carbon by land and the oceans in response to changes in atmospheric CO₂ concentration to be quantified. These fluxes can be used in conjunction with the prescribed changes in atmospheric CO₂ burden to determine the rate of fossil fuel emissions that are compatible with the prescribed CO₂ concentration pathway (see section 2a). When an ESM is forced with observed historical CO₂ concentration, the compatible fossil fuel (FF) emissions can be compared with historical anthropogenic FF emissions as a test of the model's ability to reproduce various components of the historical carbon budget realistically. For future scenarios, compatible emissions from ESMs can be used to evaluate the much less complex integrated assessment models; a good correlation between IAM emissions and ESM

emissions provides evidence that the IAM's carbon cycle-climate submodel, which generates the CO₂ concentration pathway from its emissions, is broadly consistent with the ESM's much more sophisticated carbon cycle. While an IAM provides one "realization" of future emissions compatible with a CO₂ pathway, with the different responses of their land and ocean carbon cycle components, ESMs provide a range of emissions that are compatible with a given CO₂ pathway. This range highlights the uncertainty in our understanding of the response of the coupled climate-carbon cycle system. It provides, therefore, an estimate of uncertainty in the compatible emissions rate, cumulative emissions, and the FF contribution to associated carbon budgets; for example, emissions remaining before 1.5° or 2.0°C of warming are reached.

In section 2, we discuss in more detail the way in which the compatible emissions are diagnosed, as well as the SSP scenarios. We present the compatible emissions from the nine ESMs' historical and SSP experiments in section 3 and examine which scenarios require emissions to reduce to below zero, and the timing at which this occurs. We show the ESMs' land and ocean carbon cycle responses, as well as the cumulative airborne fraction of emissions. The nine CMIP6 ESMs included in this study are summarized in Tables 1–3 and are described in further detail in section 1 of the online supplemental material.

## 2. Methods

### a. Calculation of compatible fossil emissions

Although the current generation of Earth system models typically contain a fully interactive carbon cycle in which the model can simulate the atmospheric CO₂ concentration in response to a scenario of CO₂ emissions, when taking part in model intercomparison projects, such models typically use
<table><caption>TABLE 2. Summary of CNRM-ESM2.1, IPSL-CM6A-LR, and MIROC-ES2L. See supplementary information for detailed model descriptions and explanation of acronyms used in the table.</caption><thead><tr><th></th><th>CNRM-ESM2.1</th><th>IPSL-CM6A-LR</th><th>MIROC-ES2L</th></tr></thead><tbody><tr><td>Atmosphere<br>Resolution, vertical levels</td><td>ARPEGE-Climat v6.3<br>1.4°, 91</td><td>CCSR-NIES<br>2.5° × 1.3°, 79</td><td>CCSR-NIES AGCM<br>2.81°, 40</td></tr><tr><td>Land<br>Dynamic vegetation<br>Natural vegetation PFTs</td><td>ISBA-CTRIP<br>No<br>8 trees, 3 grasses, shrub, wetlands</td><td>ORCHIDEE v2.0<br>No<br>8 trees, 4 grasses</td><td>MATSIRO + VISIT-e<br>No<br>11 vegetation types (primary and secondary)</td></tr><tr><td>Soil carbon pools</td><td>3</td><td>3</td><td>3</td></tr><tr><td>Litter carbon pools</td><td>4</td><td>4</td><td>3</td></tr><tr><td>Land-use PFTs</td><td>3 crops</td><td>2 crops</td><td>1 crop, 1 pasture</td></tr><tr><td>Net or gross transitions</td><td>Net</td><td>Net</td><td>Gross</td></tr><tr><td>Harvest</td><td>No</td><td>Wood harvest</td><td>Crop and wood harvest</td></tr><tr><td>Grazing</td><td>No</td><td>No</td><td>Yes, on pasture and rangeland</td></tr><tr><td>Fire-carbon interactions</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>Nitrogen cycle</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>Nitrogen limitation of<br>carbon uptake</td><td>Implicit N limitation: LAI reduces<br>as CO<sub>2</sub> increases, along with<br>downregulation of photosynthesis</td><td>Downregulation of<br>photosynthetic capacity</td><td>Downregulation of<br>photosynthetic capacity</td></tr><tr><td>Ocean<br>Resolution</td><td>NEMO v3.6<br>1° up to 1/3° in the tropics</td><td>NEMO v3.6<br>1° up to 1/3° in the tropics</td><td>COCO<br>1°, finer at the equator and in<br>the Arctic</td></tr><tr><td>Ocean biogeochemistry</td><td>PISCESv2-gas</td><td>PISCES-v2</td><td>OECO2</td></tr></tbody></table>

prescribed CO₂. This makes it simpler to ensure that the resulting spread of climate projections are in response to identical CO₂ forcing. However, since the atmospheric CO₂ is prescribed, and the land and ocean carbon cycle submodels

calculate how much CO₂ is taken up by land and by the oceans, it is possible to infer the rate at which additional CO₂ must be emitted to maintain the evolving atmospheric CO₂ pathway, offsetting that taken up by the land and ocean sinks.

<table><caption>TABLE 3. Summary of MPI-ESM1.2-LR, NorESM2-LM, and UKESM1-0-LL. See supplementary information for detailed model descriptions and explanation of acronyms used in the table.</caption><thead><tr><th></th><th>MPI-ESM1.2-LR</th><th>NorESM2-LM</th><th>UKESM1-0-LL</th></tr></thead><tbody><tr><td>Atmosphere<br>Resolution, vertical levels</td><td>ECHAM6.3<br>1.8° × 1.8°, 47</td><td>CAM6-Nor<br>1.9° × 2.5°, 32</td><td>UM11.2<br>1.875° × 1.25°, 85</td></tr><tr><td>Land<br>Dynamic vegetation<br>Natural vegetation PFTs</td><td>JSBACH3.2<br>Yes<br>4 trees, 2 shrubs, 2 grasses</td><td>CLM5<br>No<br>8 trees, 3 shrubs, 4 grasses</td><td>JULES-ES v1.0<br>Yes<br>5 trees, 2 shrubs, 2 grasses</td></tr><tr><td>Soil carbon pools</td><td>10</td><td>3</td><td>4</td></tr><tr><td>Litter carbon pools</td><td>8</td><td>4</td><td>None, litter flux from vegetation to soil</td></tr><tr><td>Land-use PFTs</td><td>2 crop, 2 pasture</td><td>6 crop types</td><td>2 crop, 2 pasture</td></tr><tr><td>Net or gross transitions</td><td>Gross</td><td>Net</td><td>Net</td></tr><tr><td>Harvest</td><td>Crop and wood harvest</td><td>Crop and wood harvest</td><td>Crop harvest</td></tr><tr><td>Grazing</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>Fire-carbon interactions</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Nitrogen cycle</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Nitrogen limitation of<br>carbon uptake</td><td>Yes</td><td>FUN model simulates<br>carbon cost of N uptake</td><td>NPP is downregulated</td></tr><tr><td>Ocean<br>Resolution</td><td>MPIOM1.6<br>1.5°, finer close to Antarctica<br>and Greenland, L40°</td><td>BLOM<br>1°, refined at equator</td><td>NEMO v3.6<br>1°</td></tr><tr><td>Ocean biogeochemistry</td><td>HAMOCC6</td><td>iHAMOCC</td><td>MEDUSA2.1</td></tr></tbody></table>
As described in appendix 1 of Jones et al. (2013), any increase to the combined atmosphere, land and ocean carbon store over time must be equal to the flux of “new” carbon added to the system by the combustion of fossil fuels and other processes that release to the atmosphere carbon that was previously locked away on geological time scales, such as gas flaring and cement manufacture (when using the term “fossil fuel,” we include all of these sources). By contrast, anthropogenic emissions from land-use practices release to the atmosphere CO₂ that is continually cycling through the land, ocean, and atmosphere carbon stores, on time scales from weeks to centuries, and so represent no net increase to the combined global carbon stores.

Equation (1), below, adapted from Eq. (1) of Friedlingstein et al. (2019) shows the relation between growth of the atmospheric carbon store ($G_{ATM}$), the gross land carbon sink before accounting for emissions due to land-use change ($S_{LAND}$), the ocean carbon sink ($S_{OCEAN}$), and anthropogenic emissions of CO₂ from fossil sources ($E_{FF}$), and from land-use change ($E_{LUC}$):

$$E_{FF} + E_{LUC} = G_{ATM} + S_{OCEAN} + S_{LAND}. \quad (1)$$

The budget imbalance term of Eq. (1) of Friedlingstein et al. (2019) is a residual term required to close the budget since all other terms are based on observations, and are therefore inexact. In the context of an individual model’s simulated carbon stores and the fluxes between them, all terms are known exactly, so an imbalance term is not necessary. Rearranging Eq. (1) gives

$$E_{FF} = G_{ATM} + S_{OCEAN} + (S_{LAND} - E_{LUC}). \quad (2)$$

We calculate the annual atmospheric CO₂ growth rate ($G_{ATM}$) in gigatons of carbon (1 GtC = 1 PgC = $10^{15}$ gC) per year, from the prescribed CO₂ concentration in parts per million (ppm) using the equivalence 1 part per million of CO₂ in the atmosphere is equal to 2.124 GtC (Ballantyne et al. 2015). The annual mean ocean carbon sink ($S_{OCEAN}$) is calculated from the global total of the air-to-sea flux of CO₂. The annual mean net land sink after accounting for emissions from land use ($S_{LAND} - E_{LUC}$) is equal to the global total of the net biosphere production: the gross land carbon sink ($S_{LAND}$) reduced by emissions due to land-use change, harvest, grazing, and fire ($E_{LUC}$). The net biosphere production (NBP) diagnostic, nbp, and the air-to-sea CO₂ flux diagnostic, fgco2, are two of the variables submitted to the Earth System Grid Federation (ESGF) CMIP6 data archive.<fnref n="2" /> Thus, substituting the annual mean growth rate of the prescribed CO₂ burden, fgco2, and nbp into (2) we calculate the annual mean fossil fuel emissions rate ($E_{FF}$) shown in section 3.

b. Calculation of emissions due to land-use and land-cover change

Emissions of CO₂ associated with wood harvesting and land-use and land-cover changes (LULCC) represent the greatest source of uncertainty in the global carbon cycle (Friedlingstein

et al. 2019). These emissions cannot be measured directly; observations can be made only of the net flux of CO₂ between the land and the atmosphere, that is, the natural land carbon sink minus carbon lost due to LULCC emissions. Isolating the latter and its associated uncertainty has been the focus of enormous effort for many years (Pongratz et al. 2014). The last two Global Carbon Budget estimates of LULCC emissions (Le Quéré et al. 2018; Friedlingstein et al. 2019) have been calculated as the mean of the results from two bookkeeping models (Houghton and Nassikas 2017; Hansis et al. 2015). Bookkeeping models first consider literature-derived estimates of the changes in soil and vegetation carbon stores, before and after a change in land use such as replacement of forest with crops. They then distribute the resulting carbon stock changes in time according to observation-based curves of decay of plant and soil carbon in the field or in products, and regrowth following agricultural abandonment or wood harvesting.

A multitude of terms has been used in the literature to refer to emissions from land use, with contributions from a great range of constituent processes. The scope of processes that can be simulated depends on the nature of the type of model in question, with bookkeeping models, dynamic vegetation models and fully coupled Earth system models differing in their capacity to include them (Pongratz et al. 2014). In the experiments considered here, since the atmospheric CO₂ concentration is prescribed, the CO₂ emitted by LULCC is not added to the atmosphere store.

The nine ESMs considered in the present study all include some representation of land-use change, of varying degrees of complexity, with associated CO₂ emissions. Land management processes, such as crop harvest, wood harvest, and grazing of pasture, are included in some models but not others, as summarized in Tables 1–3. The land-use emissions in a single experiment will include the contribution from the processes specific to each model, output by the appropriate combination of diagnostics; this approach to quantifying LULCC emissions is described as method S (single simulation) in Pongratz et al. (2014).

A more accurate method of determining the LULCC emissions $E_{LUC}$ is to calculate the difference in net biosphere production between a pair of simulations, one with land use changing over time, and the other with fixed land use [the experimental design denoted E2 in Pongratz et al. (2014)]. For the historical period, this is given in Eq. (3), where the two experiments are the standard CMIP6 historical, and the fixed, preindustrial land-use variant, “hist-noLu”, included in the Land Use Model Intercomparison Project (LUMIP; Lawrence et al. 2016):

$$E_{LUC} = \text{NBP}_{\text{hist-noLu}} - \text{NBP}_{\text{historical}}. \quad (3)$$

This approach identifies the total difference in land carbon store due solely to land-use change and associated processes in the historical experiment. It takes into account the additional carbon uptake as natural vegetation is allowed to regrow on abandoned agricultural land, which is neglected by the online diagnosis of land-use emissions within a single simulation. It also accounts for the loss of the additional sink capacity (LASC) that natural vegetation would have provided had it not

<footnote n="2">² https://esgf-node.llnl.gov/search/cmip6/.</footnote>
been replaced by managed land. The inclusion of the loss of additional sink capacity constitutes a key difference to estimates from the bookkeeping approach, which does not consider transient changes in environmental conditions (Pongratz et al. 2014). By using only the diagnostic of net biosphere production *nbp*, the “E2” approach does not rely on identifying the combination of diagnostics specific to each model that constitute its total LULCC emissions. However, it does require the additional computational expense of performing the no-land-use change *noLu* variant of experiments with land-use change.

### c. Scenarios

Under the CMIP5 modeling framework (Taylor et al. 2012), the response of the Earth system to realistic scenarios of future forcing was explored via the representative concentration pathways (RCPs; Vuuren et al. 2011). These consist of four distinct pathways leading to radiative forcing at 2100 of approximately 2.6, 4.5, 6.0, and 8.5 W m⁻², known as RCP2.6, RCP4.5, RCP6.0, and RCP8.5, respectively. Each pathway provides the greenhouse gas concentration and land-use scenarios, such as the area of land devoted to crop and pasture, required as forcing by Earth system models.

For CMIP6, there are eight future scenario experiments that make up ScenarioMIP (O'Neill et al. 2016); four share similar 2100 forcing levels with the RCPs, while four more lead to complementary radiative forcing levels of 1.9, 3.4 (two scenarios), and 7.0 W m⁻². The basis for each scenario is one of five SSPs, numbered 1 to 5 (Riahi et al. 2017). The SSPs describe five future narratives of economic development (Dellink et al. 2017), population change (Samir and Lutz 2017), and urbanization (Jiang and O'Neill 2017). They are designed to encompass a range of potential futures characterized by how challenging the mitigation of climate change, and adaptation to it, will be in those futures.

The scenarios are generated by intermediate-complexity integrated assessment models, which often comprise many individual submodel components (Riahi et al. 2017). They take as their inputs a wide range of variables such as scenarios of economic development, population growth, the degree of international cooperation, policy decisions on the importance of investing in education, reducing inequality, sustainability of energy resources and many other complex traits of modern society. The models are run into the future to determine how energy use and emissions of CO₂ and other greenhouse gases evolve in response to these input scenarios. For ScenarioMIP, the emissions of different greenhouse gases from the IAMs were converted into their respective concentrations using version 7.0 of MAGICC, a carbon cycle-climate model of intermediate complexity (Meinshausen et al. 2020, 2011).

The ScenarioMIP experiments are labeled SSPx-y.z where x is the SSP number from 1 to 5, and y.z is the forcing level in W m⁻², as was the case with the RCPs. The CMIP6 experiment identifier for each of the SSP experiments is sspxyz, for example, ssp370, with the exception of ssp534-over, which is an overshoot scenario, branching from ssp585 in 2040. The experiments are ordered by priority into tier 1 (ssp126, ssp245, ssp370, ssp585) and tier 2 (ssp119, ssp434, ssp460, ssp534-over).

We present results from all eight SSP experiments, whose CO₂ concentration pathways are shown in Fig. 1a, with the global population in Fig. 1b. The crop and pasture fractions of the land-use scenarios are shown in Fig. 1c, with forest area in Fig. 1d.

### d. Earth system models

The nine Earth system models included in the present study, with the abbreviated form by which they will be referred in the text in parentheses where necessary, are ACCESS-ESM1.5 (ACCESS), CanESM5, CESM2, CNRM-ESM2.1 (CNRM), IPSL-CM6A-LR (IPSL), MIROC-ES2L (MIROC), MPI-ESM1.2-LR (MPI), NorESM2-LM (NorESM), and UKESM1-0-LL (UKESM). The models are described in detail in supplementary section 1 and summarized in Tables 1–3. All nine ESMs ran the historical experiment, which is a requisite for all models participating in CMIP6, and the four tier 1 SSP experiments. The tier 2 SSP experiments were performed by 3, 4, or 5 ESMs, as summarized in Table 4.

### e. Reference data

Where possible, the historical plots and tables include reference data from the Global Carbon Budget, 2019 (GCB2019; Friedlingstein et al. 2019). The global fossil emissions rate (from fossil carbonates as well as fossil fuels) $E_{FF}$ are from Gilfillan et al. (2019). The $E_{LUC}$ emissions from GCB2019 are the average of updated versions of the data from Houghton and Nassikas (2017) and Hansis et al. (2015). The terrestrial and ocean carbon sinks, $S_{LAND}$ and $S_{OCEAN}$, respectively, come directly from the GCB2019. The $E_{FF}$ is based on historical records of fossil emissions, but the $E_{LUC}$, $S_{OCEAN}$, and $S_{LAND}$ are all derived from models constrained by observations where possible.

## 3. Results and discussion

The results have been divided into a historical section from 1850 to 2014 (section 3a) and a future section covering all SSP experiments from 2015 to 2100 (section 3b). The time series shown were constructed from ensemble means, where possible, of the *nbp*, *fgco2*, *cVeg*, *cSoil*, *cLitter* (vegetation, soil, and litter carbon store), *fharvest* (harvest flux), and *tas* (surface air temperature) diagnostics submitted to ESGF. The models' individual responses to the SSPs for each variable are shown in the online supplemental material, which also contains tables summarizing ocean and land carbon uptake historically and under the SSPs.

### a. Historical

#### 1) COMPATIBLE FOSSIL FUEL EMISSIONS

Equation (2) has been used to calculate time series of compatible fossil fuel CO₂ emissions for the ESMs' historical experiment (Fig. 2). For comparison, the data-based estimates of historical emissions from fossil fuel burning and cement production (Gilfillan et al. 2019) are shown (black dashed lines). The multimodel mean follows very closely the observed historical emission rate, except for a period of around 15 years,
<figure><img src="image_2.png" /><figcaption>FIG. 1. Scenarios of (a) atmospheric CO₂ concentration (parts per million), (b) global population (billions of people), and area (millions of hectares) of the land surface covered by (c) crop and pasture and (d) forest, of the eight SSPs. Historical CO₂ concentration is shown in the inset in (a).</figcaption></figure>

from the early 1940s. After decades of rising, in 1941 CO₂ began to stabilize and remained between 310 and 312 ppm from 1941 to 1951 before rising again [Etheridge et al. (1996), confirmed more recently by MacFarling Meure et al. (2006) and Rubino et al. (2013)]. Hence, during this period of the historical experiment, G<sub>ATM</sub> reduced to almost zero, leading the ESMs to simulate a reduction in E<sub>FF</sub>, in the absence of any compensating factor to increase the ocean or land sinks. In reality, the data-based fossil fuel emissions rate showed no such lull during the observed plateau in concentration at this time,

increasing from 1.33 to 1.76 GtC yr⁻¹ from 1941 to 1951. This implies the existence in the real world of an enhanced carbon sink during this period that has not yet been explained fully by the observational datasets available (Bastos et al. 2016). The ocean is likely to have played a role, with enhanced uptake by the Southern and tropical Pacific Oceans associated with the strong El Niño of the early 1940s (Joos et al. 1999). There could have been additional uptake on land due to land-use processes and changes not accounted for in the LUC datasets, nor therefore in estimates of land uptake of the time. Due to the

<table><caption>TABLE 4. Summary of the ScenarioMIP SSP experiments performed by each model.</caption><thead><tr><th></th><th colspan="7">SSP</th></tr><tr><th></th><th>1-1.9</th><th>1-2.6</th><th>2-4.5</th><th>3-7.0</th><th>4-3.4</th><th>4-6.0</th><th>5-3.4over</th><th>5-8.5</th></tr></thead><tbody><tr><td>ACCESS-ESM1.5</td><td>—</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>—</td><td>—</td><td>✓</td></tr><tr><td>CanESM5</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>CESM2</td><td>—</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>—</td><td>—</td><td>✓</td></tr><tr><td>CNRM-ESM2.1</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>IPSL-CM6A-LR</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>MIROC-ES2L</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>—</td><td>—</td><td>✓</td></tr><tr><td>MPI-ESM1.2-LR</td><td>—</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>—</td><td>—</td><td>✓</td></tr><tr><td>NorESM2-LM</td><td>—</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>—</td><td>—</td><td>✓</td></tr><tr><td>UKESM1-0-LL</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>✓</td><td>✓</td></tr><tr><td>Total</td><td>5</td><td>9</td><td>9</td><td>9</td><td>4</td><td>3</td><td>4</td><td>9</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Time series of compatible fossil fuel emissions (GtC yr⁻¹) from 1850 to 2014 simulated by the nine ESMs. Solid black line is the model mean. The dashed black and white line is the data-based fossil fuel rate from the Gilfillian et al. (2019) as cited in the GCB2019 (Friedlingstein et al. 2019), with decadal mean emission rates from the 1960s to the 2010s and associated uncertainty, indicated by the white error bars, from GCB2019. (b) Cumulative fossil fuel emissions (GtC) for the same period.</figcaption></figure>

Second World War, the 1940s were characterized by enormous socioeconomic upheaval with changes to land management, food production, and relocation of population; not all changes contributing to LUC emissions were recorded accurately, such as abandonment of land and changes to wood harvest (Bastos et al. 2016). It should also be noted that this era predates direct measurement of atmospheric CO₂ concentration, and estimates of fossil fuel emissions were less accurate than they are today. So our understanding of the 1940s CO₂ plateau and of the fluxes potentially contributing to the implied sink remains incomplete (Bastos et al. 2016). It would be wrong to conclude, therefore, that because the ESM's fossil emissions fall short of the observations-based emissions during this period, the models systematically underestimate FF emissions during periods of constant or very slowly changing CO₂, as occur in some of the SSP scenarios discussed in section 3b.

The ESMs' temperature response to the historical forcing is shown in Fig. S1 in the online supplemental material.

## 2) CUMULATIVE FOSSIL FUEL EMISSIONS

Figure 2b shows the cumulative fossil fuel emissions from the ESMs, with the data-based estimate shown for comparison (Gilfillan et al. 2019). The mean and standard deviation of the total FF emissions from 1850 to 2014 was 394 ± 59 GtC, compared to observations-based estimate of 400 ± 20 GtC for the same period (Gilfillan et al. 2019). The model mean is therefore very close to the observed estimate, lying well within its uncertainty range (Table 5); however, exclusion of the highest two model values results in a model mean of 367 GtC, below the lower limit (380 GtC) of the observed range. The ESMs' underestimate of emissions during the period of CO₂ plateau is largely responsible for this; cumulative totals of

<table><caption>TABLE 5. Cumulative compatible fossil fuel CO₂ emissions (GtC) simulated by the ESMs for the historical period (1850–2014) and for all SSPs (from 1850 to 2100, including the historical emissions). The Obs/IAM row contains the cumulative emissions from 1850 to 2100, incorporating the observed historical emissions from 1850 to 2014 (Gilfillan et al. 2019); uncertainty range (±SD) in the historical figure comes from Friedlingstein et al. (2019)</caption><thead><tr><th rowspan="2"></th><th colspan="9">SSP</th></tr><tr><th>1850–2014</th><th>1-1.9</th><th>1-2.6</th><th>2-4.5</th><th>3-7.0</th><th>4-3.4</th><th>4-6.0</th><th>5-3.4over</th><th>5-8.5</th></tr></thead><tbody><tr><td>ACCESS-ESM1.5</td><td>489</td><td>—</td><td>814</td><td>1238</td><td>1841</td><td>—</td><td>—</td><td>—</td><td>2515</td></tr><tr><td>CanESM5</td><td>349</td><td>526</td><td>766</td><td>1309</td><td>2054</td><td>706</td><td>1493</td><td>911</td><td>2875</td></tr><tr><td>CESM2</td><td>339</td><td>—</td><td>794</td><td>1240</td><td>1870</td><td>—</td><td>—</td><td>—</td><td>2609</td></tr><tr><td>CNRM-ESM2.1</td><td>506</td><td>735</td><td>950</td><td>1450</td><td>2083</td><td>978</td><td>1598</td><td>1188</td><td>2798</td></tr><tr><td>IPSL-CM6A-LR</td><td>393</td><td>614</td><td>822</td><td>1261</td><td>1878</td><td>797</td><td>1403</td><td>987</td><td>2538</td></tr><tr><td>MIROC-ES2L</td><td>406</td><td>581</td><td>800</td><td>1325</td><td>2085</td><td>—</td><td>—</td><td>—</td><td>2745</td></tr><tr><td>MPI-ESM1.2-LR</td><td>352</td><td>—</td><td>782</td><td>1226</td><td>1885</td><td>—</td><td>—</td><td>—</td><td>2574</td></tr><tr><td>NorESM2-LM</td><td>346</td><td>—</td><td>817</td><td>1266</td><td>1894</td><td>—</td><td>—</td><td>—</td><td>2646</td></tr><tr><td>UKESM1.0-LL</td><td>364</td><td>507</td><td>701</td><td>1145</td><td>1753</td><td>673</td><td>—</td><td>844</td><td>2432</td></tr><tr><td>Mean ± SD</td><td>394 ± 59</td><td>593 ± 81</td><td>805 ± 62</td><td>1273 ± 79</td><td>1927 ± 111</td><td>789 ± 118</td><td>1498 ± 80</td><td>983 ± 129</td><td>2637 ± 136</td></tr><tr><td>Obs/IAM</td><td>400 ± 20</td><td>524</td><td>726</td><td>1212</td><td>1909</td><td>779</td><td>1366</td><td>906</td><td>2580</td></tr></tbody></table>
several of the ESMs diverge from the observations-based total at this time and remain so until 2014 (Fig. 2b). Only two of the nine ESMs (IPSL and MIROC) lie within the observed cumulative 1850–2014 range (Table 5).

Before looking at how the land and ocean sinks contribute to the historical FF emission rate for the nine ESMs in detail, it is worth examining a subset of simulations with fixed land use, to see the extent to which land-use change influences the magnitude of, and model spread in, compatible fossil fuel emissions.

### 3) INFLUENCE OF HISTORICAL LAND-USE EMISSIONS ON FOSSIL FUEL EMISSIONS

Of the nine ESMs included in this study, seven have submitted output from both *historical* and *hist-noLu* experiments to the CMIP6 archive. The LULCC emissions rate, calculated using Eq. (3), is shown in Fig. 3c, with the cumulative total in Fig. 3d, for comparison with the equivalent FF emissions (Figs. 3a,b, respectively) from the *historical*, with GCB2019 reference values for comparison. Note that the GCB2019 LULCC emission estimate is based on the bookkeeping approach, which does not include loss of additional sink capacity. The LASC adds approximately another 20% to the LULCC emissions estimate [GCB2019 states 10%, but this was revised in Gasser et al. (2020)]. For clarity, larger versions of the individual panels of Fig. 3 appear in the supplemental material as Figs. S13–S15.

The spread in the magnitude of the LULCC emissions among the seven models is immediately apparent, with CESM2's, which are very close to the observed estimates, being more than eight times those from CNRM. However, it is noticeable that all ESMs exhibit lower cumulative $E_{LUC}$ than GCB2019, with the exception of the last few years of MPI, despite their accounting for the loss of additional sink capacity. This hints toward LUC processes missing from the ESMs that are included in the GCB2019 estimate due to the construction of the bookkeeping models, which use observed carbon densities for various managed and unmanaged vegetation types.

The model spread in LULCC emissions is mirrored approximately by the spread in FF emissions, cumulatively (Figs. 3d,b, respectively); greater LULCC emissions are paired with lesser FF emissions, as is seen most clearly with those ESMs at the edges of the model spread, CNRM, MPI, and CESM2. Broadly speaking, the greater the extent to which a model's land carbon uptake is limited by land use, the lesser the amount of fossil CO₂ can be emitted for a given CO₂ concentration pathway. Combining the FF and LULCC emissions (Fig. 3e) significantly reduces the model spread (Fig. 3f) as the influence of land use on both LULCC and FF emissions cancels out, so differences only in the models' natural carbon sink are in evidence; the model spread at 2013 for FF + LULCC is 77 GtC, compared to 168 GtC for FF and 173 GtC for LULCC. The differences in the models' natural carbon sink are therefore dwarfed by the differences between the net land carbon sink after accounting for land-use practices (as will be seen in Fig. 5b). For the seven models included in Fig. 3, the presence of land use is therefore the greatest source of uncertainty in future FF emissions, and on future carbon budgets utilizing these.

Similar analysis cannot be done for the SSP experiments since the *noLu* variant of the experiments have not been

requested under CMIP6 protocols for ScenarioMIP, C4MIP or LUMIP.

The magnitude of LULCC emissions calculated from the E2 paired experiments depends on several factors, including the spatial extent and dynamics of the land-use forcing dataset used in the *historical*; the carbon density of the vegetation being cleared; the land-use practices such as irrigation, fertilization, harvesting, and grazing, which impact the agricultural vegetation; and the responses of both the natural and agricultural vegetation to the changing climate and CO₂ concentration. The response of natural vegetation to future climate and CO₂ will depend on whether the vegetation distribution is fixed, or a dynamic vegetation model is utilized that allows the vegetation cover to migrate to latitudes and regions that become habitable as the climate changes. Land management practices that transfer living biomass to product pools, such as wood harvest and crop harvest, also play a role. For the historical experiment, although the spatial land-use forcing is provided in each case by the Land Use Harmonization 2 (LUH2) dataset (Hurtt et al. 2020), the combination of land-use states from LUH2 used is model dependent. So too is the method by which vegetation types within a grid cell change. “Net transitions,” considering only the net differences of each vegetation type’s fractional coverage following a land-cover change, can underestimate the impact on carbon stores compared to the “gross transitions” approach. Under gross transitions, the localized carbon impact at subgrid scale of each individual transition between neighboring vegetation types is considered explicitly; for example, primary forest of a young age structure replaces abandoned cropland, which is distinct from older, secondary forest already present elsewhere in the gridbox (Yue et al. 2018). Of the nine ESMs in the present study, two (MIROC, MPI) consider gross transitions; these, along with CESM2, get closest to the GCB2019 cumulative $E_{LUC}$.

### 4) LAND CARBON RESPONSE

The simulated historical *nbp*, the land carbon uptake after accounting for carbon lost due to land-use change, is shown in Fig. 4. The GCB2019 provides observationally constrained multimodel estimates of the land carbon sink in each of the decades from the 1960s onward; these are included in Fig. 4 and are compared to the model's means for the same time periods in Table S2. The model mean sits well within the uncertainty range of the decadal mean estimates. The GCB2019s $S_{LAND}$ and $E_{LUC}$ estimates of 185 ± 50 GtC and 195 ± 75 GtC, respectively, result in net land carbon uptake from 1850 to 2014 of -10 ± 125 GtC. The ESMs' net uptake ranges from -35 to +142 GtC over this period, with a mean of +20 GtC, as summarized in Table S3.

The nine ESMs' preindustrial (1850) land carbon stores, with partitioning between soil, litter, and vegetation, are shown in Fig. 5a. The change in land carbon (cumulative *nbp*) relative to 1850 is shown in Fig. 5b, with changes in vegetation carbon (*cVeg*) and soil carbon (*cSoil* + *cLitter*) in Figs. 5c and 5d, respectively. Broadly, the models' historical land carbon uptake is dominated by the change in vegetation carbon, with the initial reduction due to deforestation offset in later decades by
<figure><img src="image_4.png" /><figcaption>FIG. 3. (left) Time series and (right) cumulative totals of CO₂ emissions from seven ESMs with observational fossil fuel emissions from Gilfillan et al. (2019) and land-use emissions from the GCB2019, which also provides the uncertainty range in cumulative totals at 2014 (Friedlingstein et al. 2019). (a) Compatible fossil fuel emissions (GtC yr⁻¹) calculated for the *historical* experiment. (b) Cumulative fossil fuel emissions from the *historical* experiment. (c) Land-use emissions calculated as land carbon uptake (cumulative NBP) in *hist-noLu* minus that in *historical*. (d) Cumulative land-use emissions. (e) Combined fossil fuel and land-use emissions. (f) Cumulative total of combined fossil fuel and land-use emissions. Note that (c) to (f) contain two data-based estimates of E_LUC. The thick black curve is the dataset of land-use emissions provided under the C4MIP protocols for driving the emissions-driven historical *esm-hist* experiment for models that do not have the ability to simulate LUC emissions interactively with the atmosphere (Jones 2016). These are derived from Hansis et al. (2015) and Houghton et al. (2012). The thin black curve is E_LUC from Friedlingstein et al. (2019), the cumulative total calculated using Table 8 therein.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time series of net biosphere production (GtC yr⁻¹) from 1850 to 2014 calculated by the ESMs. Data-based estimates of decadal mean net land carbon uptake for the 1960s, 1970s, 1980s, 1990s, and 2000s and from 2009 to 2018 are overlaid as squares with uncertainty range, from Table 6 of Friedlingstein et al (2019).</figcaption></figure>

the increase due to enhanced uptake by CO₂ fertilization and agricultural use of nitrogenous fertilizer. ACCESS is the only model with continued increase in both vegetation and soil carbon. Although ACCESS does include land-use change, its vegetation carbon is approximately constant until the 1960s when it begins to rise in common with the other ESMs. This is because crops mostly replace natural grasses, with some cleared forest in the tropics balanced by growth of new forest at northern mid- to high latitudes, resulting in little change in vegetation carbon overall until the latter half of the century.

The multimodel mean change in soil carbon shows a very gradual increase throughout the historical period. By contrast, MIROC's soil carbon reduces until the final few decades, when it begins to increase. In MIROC, replacement of primary/ secondary vegetation by crop and pasture of lower productivity reduces the soil and litter carbon until they are balanced by the lower litter inputs. In the latter half of the twentieth century, due to the effects of CO₂ fertilization and nitrogen fertilizer application, vegetation, litter, and soil carbon increase. Considering MIROC's intermediate CO₂ concentration-carbon feedback strength (Arora et al. 2020), the relatively steep increase of land carbon in this period is likely induced by accumulation in cropland, as suggested in Hajima et al. (2020).

The presence of a crop harvest flux can impact the vegetation carbon and/or soil carbon stores significantly. The crop harvest from five ESMs' historical simulations relative to their piControl is shown in Fig. 6, as a rate and cumulatively. The cumulative total at 2014 ranges from 60 to 246 GtC, so represents a fairly large impact on vegetation, litter, or soil carbon stores depending on how harvest is configured in each model. In UKESM, for example, 30% of the litter flux from crops to soil is intercepted to represent harvest, so UKESM's increase in historical soil carbon is not as great as it would otherwise be, but its living vegetation carbon store is not directly affected by the harvest.

Of the three ESMs that do not include a crop harvest, two (ACCESS, CNRM) exhibit the two greatest increases in historical soil and land carbon, contributing to the two greatest fossil fuel emission rates of all nine ESMs. Conversely, the third ESM in which crop harvest is not represented, CanESM5, sees a downward trend in soil carbon throughout the historical period that steepens in the final few decades. This is because in CanESM5 the soil decomposition rate over croplands is higher and the fraction of humified litter that is transferred to the soil carbon, as opposed to respired to the atmosphere, is lower than over natural vegetation. Consequently, as natural vegetation is replaced by croplands over the historical period, a decrease in global soil carbon is obtained, as is also seen in empirical measurements (Wei et al. 2014).

During harvest in CESM2 and NorESM2 (in both of which the land component is CLM5), even though crop grain carbon is removed, the increased productivity of crops due to fertilization, irrigation, and higher productivity plants actually leads to increases in soil carbon where crops are grown. This response is the opposite to what is observed in the real world, which is possibly due to lack of a representation of tillage in the model; tillage would likely lead to increased respiration, resulting in soil carbon losses. Further, what happens with soil carbon will depend on how much plant matter is removed during harvest; CLM5 assumes only grain is removed, but in reality, usually more than just grain is removed during the harvest process (Lombardozzi et al. 2020).

CESM2 and NorESM2 also simulate a wood harvest flux, which can result in a greater loss of vegetation carbon than deforestation. The wood harvest flux has a greater effect on vegetation carbon store as the regrowth of the harvested trees occurs on a much longer time scale than is the case with crops. This also impacts the soil carbon evolution in CESM2.

MIROC and MPI include an additional carbon loss due to grazing of pasture, which in MIROC results in a reduction in leaf area and reduces productivity.

### 5) OCEAN CARBON RESPONSE

The air-to-sea CO₂ fluxes are very similar in all nine ESMs (Fig. 7a) throughout the historical period. The model mean lies slightly lower than the six decadal mean values from GCB2019, but well within their uncertainty range, as is true for all individual models; these are compared with the GCB2019 reference values in Table S4. The models sit at (ACCESS) or below the GCB2019 change in ocean carbon store of 150 GtC from 1850 to 2014 (Table S4). IPSL and CNRM fall just below the lower end of the uncertainty range of the reference values from GCB2019.

### 6) AIRBORNE FRACTION

The rate of growth of atmospheric CO₂ concentration is dependent on the airborne fraction (AF) of emitted CO₂, the proportion of emissions that remains in the atmosphere rather than being absorbed by the oceans or land. The AF varies with the rate of emissions; if uptake by the sinks is outpaced by emissions, the emitted CO₂ will begin to accumulate in the atmosphere, resulting in a higher AF and greater radiative forcing. Whether the airborne fraction of emissions has begun
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Initial land carbon stores in all ESMs at 1850, partitioned into vegetation, litter, and soil carbon. (b) Change in land carbon from 1850 to 2014 calculated as the cumulative *nbp* diagnostic. (c) Change in vegetation carbon and (d) change in +litter carbon over the same period. The values in (c) and (d) are calculated from the *cVeg* and *cSoil* + *cLitter* diagnostics, respectively.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Harvest carbon flux from the historical simulation of five ESMs. The data are annual means calculated from one ensemble member, with the average *fHarvest* from at least a 10-yr section of the *piControl* subtracted. (b) Cumulative totals of harvest fluxes shown in (a).</figcaption></figure>

to increase in recent history and whether it will do so over the coming century are therefore important questions.

Several studies in recent years have attempted to find evidence for a trend in AF from observations taken over the last 50 years or so. The airborne fraction is typically determined as a fraction of total anthropogenic emissions, from fossil fuel as well as land-use change (Canadell et al. 2007; Denman et al. 2007; Le Quéré et al. 2009; Raupach et al. 2008; Rayner et al. 2010; Bennedsen et al. 2019), though it has also been estimated from fossil fuels and other industrial sources alone, for example, from observations by Keeling et al. (1995) and for the CMIP5 models by Jones et al. (2013) and Ciais et al. (2013). These studies have suggested a weak increase in the airborne fraction of emitted carbon since the 1960s, although the robustness of this finding has been questioned (Knorr 2009; Raupach et al. 2014; Ballantyne et al. 2015). Since the focus of

this paper is on compatible fossil fuel emissions, and due to the range of approaches to land-use emissions in the CMIP6 ESMs, we calculate the AF from the fossil fuel emissions alone, and the cumulative airborne fraction (CAF) rather than the instantaneous annual mean AF. The annual mean AF is a quotient constructed from three potentially small fluxes (rise in atmospheric CO₂, *nbp*, and *fgco2* for a particular year) and can vary widely from year to year, whereas the cumulative AF is calculated as the total amount by which the atmospheric burden has increased since the preindustrial, expressed as a fraction of the cumulative fossil fuel emissions. The cumulative AF has a much stronger signal-to-noise ratio after several decades of emissions and will change over time to reflect the ability of the land and ocean sink to keep pace with the emissions.

The time series of cumulative airborne fraction for the ESMs in the historical period is shown in Fig. 8. There is little

<figure><img src="image_8.png" /><figcaption>FIG. 7. (a) Time series of air-to-sea flux of CO₂ from the nine ESM's *historical* experiments. Data-based estimates of decadal mean fluxes for the 1960s, 1970s, 1980s, 1990s, and 2000s and from 2009 to 2018 are overlaid as squares with an uncertainty range. (b) Change in ocean carbon store calculated as the cumulative total of the air-to-sea flux shown in (a). The square at 2014 indicates the data-based estimate of cumulative carbon uptake from 1850 to 2014. All observations are from Table 6 of Friedlingstein et al. (2019).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8. Time series of airborne fraction of fossil fuel emissions from the 9 ESM's historical simulation. Black squares and uncertainty ranges are observations-based estimates from Table 6.1 of Ciais et al. (2013).</figcaption></figure>

agreement between models in the trend and magnitude of the CAF historically. Four models' FF emissions are exceeded by the increment in atmospheric CO₂ resulting in a CAF exceeding unity, until the final few decades of the twentieth century; as a result, the same is true of the model mean CAF until 1940. By the end of the twentieth century, all models are converging toward a CAF of between 0.6 and 0.7, with the exception of the two ESMs in which the land carbon store increases throughout the historical period (ACCESS and CNRM). The IPCC Fifth Assessment Report gives estimates of the airborne fraction of fossil fuel emissions for the 1980s, 1990s, and two decadal periods of the present century, which also exhibit a reducing trend, from 0.62 in the 1980s to 0.51 in the first decade of this century (black boxes with uncertainty range on Fig. 8).

### 7) POTENTIAL ROLE OF NON-CO₂ GHGs AND AEROSOLS ON LAND AND OCEAN CARBON UPTAKE

Carbon uptake by the biosphere can be influenced directly by forcings other than CO₂ (Gasser et al. 2017). Aerosol-induced cooling can impact land carbon uptake directly. Zhang et al. (2019) showed that aerosol cooling increased historical land carbon uptake by 11.6 to 41.8 GtC, due primarily to increases in net biosphere production in the mid- to high latitudes; at lower latitudes, aerosol cooling reduced gross primary production (GPP) by more than total ecosystem respiration, reducing NBP.

Complex aerosol-cloud interactions are expected to lead to changes in precipitation, with implications for land carbon uptake (e.g., Fig. S11 of Zhang et al. 2019). Unger et al. (2017) used a coupled aerosol-carbon cycle-climate model to simulate the period 1996 to 2005 with and without historical aerosol levels and found that the inclusion of aerosols led to a global precipitation reduction of 0.08 mm day⁻¹ with regional increases and decreases in NPP.

A further influence of aerosols on land carbon uptake is via their increased scattering of light, increasing the proportion of diffuse to direct light and allowing it to penetrate further into the canopy, enhancing photosynthetic capacity and increasing GPP (Mercado et al. 2009).

ESMs differ markedly in their representation of aerosols and atmospheric chemistry of aerosol precursor species (e.g.,

Thornhill et al. 2021). It is beyond the scope of this paper to examine in depth these differences and the resulting aerosol concentrations, but the effective radiative forcing (ERF)<fnref n="3" /> from 1850 to 2014 due to the presence of aerosols (ERF_aer) has been calculated by Smith et al. (2020) for 6 ESMs in the present study. The ERF_aer ranges from −0.63 W m⁻² (IPSL), to −1.37 W m⁻² (CESM2), a span of 0.74 W m⁻². This illustrates that aerosol-induced climate cooling and associated changes to precipitation could vary quite strongly between models, with the land carbon uptake, and therefore the compatible fossil emissions, influenced accordingly.

The AerChemMIP (Collins et al. 2017) experimental design allows calculation of the ERF at 2014 due to well-mixed greenhouse gases (CO₂, CH₄, N₂O, and halocarbons, ERF_wmghg), CO₂ only (ERF_CO₂), and, therefore, well-mixed GHGs excluding CO₂ (ERF_nonCO₂ghg). The range of ERF_nonCO₂ghg is smaller than ERF_aer, from 0.71 to 1.15 W m⁻², a span of 0.44 W m⁻²; this is smaller than the spread in ERF_aer but shows that differences in the models' responses to the non-CO₂ greenhouse gas forcing also plays a nontrivial role in the evolution of the ESMs' climate and therefore of the carbon cycle and compatible emissions. Ozone can also affect land carbon sinks directly, through its impact on leaf health (Sitch et al. 2007).

Aerosol-induced cooling also increases solubility of gaseous CO₂ in seawater, and to some extent alters ocean circulation and biological production, all of which potentially impact ocean carbon sinks, albeit to a lesser extent than that expected

<footnote n="3">³ ERF is the radiative forcing (RF) after allowing for the adjustment of all tropospheric and land surface properties on a range of time scales (though typically referred to as “rapid adjustments”), including changes to vegetation and land ice cover, and is evaluated after these have regained equilibrium following a perturbation in forcing such as a change in CO₂ or aerosol. By contrast, RF is defined as the change in downward radiative flux at the tropopause after allowing only stratospheric temperatures to adjust to radiative equilibrium, with all other state variables held constant at the unperturbed state. The ERF represents a better indicator than RF of the eventual global mean temperature in response to a perturbation in forcing. For aerosols, in particular, ERF can be quite different from RF (Smith et al. 2020).</footnote>
<figure><img src="image_10.png" /><figcaption>FIG. 9. (a) Time series of compatible fossil fuel emissions (GtC yr⁻¹) for the eight SSP scenarios simulated by the nine ESMs. Solid lines are the model mean, the dashed lines are the fossil fuel rate generated by the IAMs that provided the CO₂ concentration pathways used to drive the ESMs. Shaded areas indicate model spread for each scenario. (b) Cumulative totals of fossil fuel emissions (GtC) shown in (a).</figcaption></figure>

from the land biosphere (Tjiputra et al. 2016; Lauvset et al. 2017). Several ESMs (MIROC, NorESM2, and UKESM) further simulate oceanic dimethyl sulfide (DMS) emissions, which is the largest natural source of atmospheric sulfur and contributes to aerosol formation and can therefore feed back to the oceanic carbon sink.

b. *Future*

1) COMPATIBLE FOSSIL FUEL EMISSIONS

For the SSP experiments, the ESM FF emissions are com-
pared with the emission rates from the IAMs in Fig. 9a. There
is generally strong agreement between the IAM and the ESMs;
the IAM FF emission rates lie within the spread of the
ESMs, with the exception of SSP4-6.0, which falls slightly
below the ESMs' range. The SSP4-6.0 was performed by
fewer ESMs (three) than any other SSP experiment; these
include CNRM and CanESM, which simulate two of the
three largest cumulative FF emissions for SSP2-4.5, the
scenario with emissions closest to those of SSP4-6.0. Had
more ESMs performed SSP4-6.0, the range of ESM emis-
sions might well have encompassed those of the IAM. (The
fossil fuel emission rates are shown for the models individ-
ually in Fig. S3.)

Three of the SSPs require the emissions to reduce to below
zero by 2100 to adhere to the CO2 pathway, in at least one
ESM. In all ESMs, SSP1-2.6 emissions become negative,
though only very slightly in NorESM2, and the same is true of
the five models that performed SSP1-1.9 and SSP4-3.4. The
interannual variability of compatible emissions varies a lot
between ESMs; some have very little (CanESM5, CNRM,
CESM2), others more so (ACCESS, MIROC). MIROC's
emissions are characterized by a fairly cyclical variability with
a period of about 10 years and an amplitude of 5 GtC yr-1. This
is likely due to the unrealistically large El Niño-Southern
Oscillation amplitude exhibited by MIROC that causes similar
variability in global temperature (Hajima et al. 2020) as can be
seen in Fig. S2.

Of interest to policymakers, particularly in the case of
overshoot scenarios, is the point at which emission rates peak
and start to decline. Figure 10 shows the peak emissions rate
versus the year in which this occurs for all nine ESMs and the
IAMs. For the lower emissions scenarios, the IAM and ESMs
are in close agreement for the magnitude of maximum emis-
sions, but the timing is more uncertain, with a range of 18 years
for SSP1-2.6 (2015-32) and 15 years for SSP2-4.5 (2043-57);
this is due in part to the broader peak and decline of the CO2
curve in the lower concentration scenarios. For the higher
emissions scenarios, SSP5-8.5 and SSP3-7.0 particularly, there
is a greater spread in the maximum emissions rate; for SSP3-7.0
emissions climb right through to 2100, so the year of highest
emissions occurs uniformly in the last few years of the simu-
lation, whereas in SSP5-8.5 emissions peak and decline, so
there is inevitably uncertainty in the year in which the peak
occurs, from 2086 to 2100.

2) CUMULATIVE FOSSIL FUEL EMISSIONS

The cumulative ESM FF emissions are compared with the emission rates from the IAMs (dashed colored lines) in Fig. 9b, with their cumulative totals from 1850 to 2100 in Table 5. The ESMs' cumulative fossil fuel emission rates are depicted individually in Fig. S4. The IAMs' cumulative totals were constructed with the FF emissions data from 1850 to 2014 (Gilfillan et al. 2019).

For the two higher emission scenarios, SSP5-8.5 and SSP3-
7.0, although the IAMs' FF emission rates increasingly drop
below the ESMs' mean toward the end of the century, periods
of higher emissions earlier mean that their cumulative totals
are in very close agreement. The IAMs also slightly underes-
timate the FF rates of the lower emission scenarios compared
to the ESMs, but the IAM cumulative total lies within the ESM
spread for all scenarios, except SSP4-6.0 for the reasons de-
scribed above. The ESM mean cumulative emissions for SSP4-
3.4 match almost exactly those constructed from historical
data-based values and the IAM and are almost identical to the
those from SSP1-2.6, although all nine ESMs performed the
<figure><img src="image_11.png" /><figcaption>FIG. 10. Maximum rate of FF emissions (GtC yr⁻¹) experienced during each SSP by each ESM vs year in which this occurred. Symbols denote SSP and ESMs are distinguished by color, with IAM in black.</figcaption></figure>

latter while only four performed SSP4-3.4. The atmospheric CO₂ is slightly higher in SSP4-3.4 than SSP1-2.6, but the former has much higher crop and pasture area (Fig. 1c), so the impact of land use brings its compatible emissions down in line with those of SSP1-2.6. Generally, where the ESMs disagree with the IAMs, the cumulative emissions from ESMs are higher, although UKESM is the exception; for SSP3-7.0 and SSP5-8.5 UKESM's cumulative emissions are lower than the emissions from IAMs by a small amount.

### 3) TIMING OF NET-ZERO EMISSIONS REQUIRED FOR 1.5° AND 2.0°C OF WARMING

SSP1-1.9 has lower radiative forcing than all of the RCPs of CMIP5 and was introduced as a result of the Paris Agreement of 2015 to explore the potential for limiting warming to 1.5°C above the preindustrial level (Rogelj et al. 2018). Based on RCP2.6, SSP1-2.6 was designed to be consistent with warming of the order of 2.0°C or lower. The assumptions underlying both of these low warming scenarios include extensive

investment in sustainability and a reduction in fossil fuel use, leading to net zero fossil fuel emissions during the second half of the twenty-first century, and negative emissions subsequently. Once such an ambitious warming limit has been decided upon, a further decision to make is the target year by which emissions should reduce to zero to achieve this. We will now examine the time frame within which the models predict fossil fuel emissions must drop to zero to adhere to the CO₂ pathway, and the change in global mean temperature the models exhibit under each.

Figure 9a shows that in all ESMs, the emission rate drops to below zero for SSP1-1.9, SSP1-2.6, and SSP4-3.4, the three SSPs with the lowest CO₂ forcing, as well as SSP5-3.4over, the overshoot scenario branching from SSP5-8.5 at 2040 with intensive mitigation. Table 6 shows the year in which the fossil fuel emissions rate first reduced to zero, when smoothed by an 11-yr running mean.

For SSP1-1.9, emissions first become zero in 2056 for the IAM; the five ESMs that ran that scenario all require emissions to become negative, starting in 2056 (UKESM) through to 2071 (CNRM). In the case of SSP1-2.6, all nine ESMs see emissions becoming negative during the period 2076–86, compared to 2076 for the IAM. Only four ESMs ran SSP4-3.4, all requiring a reduction to zero from 2080 to 2091, enclosing the IAM's value of 2084. The overshoot scenario, SSP5-3.4over, sees emissions in the ESMs becoming zero in 2068, as is true of the IAM, through to 2078.

The temperature response of the ESMs to the SSPs is shown as a model mean in Fig. 11, and individually in Fig. S2. The multimodel-mean warming at 2100 relative to the 1850–99 average ranges from 1.5° to 2.2°C for the five models that performed SSP1-1.9, with a mean of 1.9°C (Table S1). The mean of the nine ESMs' warming at 2100 under SSP1-2.6 is 2.0°C, with a range of 1.1°C (NorESM2) to 2.8°C (CanESM5 and UKESM); five ESMs exhibit warming at 2100 of 2.0°C or less under this scenario.

### 4) LAND CARBON RESPONSE

The ESMs' land carbon response to the scenarios is the net result of a combination of their natural carbon cycle behavior as the CO₂ and climate change, and the extent to which this is perturbed by the prescribed land-use forcing, and the land management practices implemented within the ESM. In the

<table><caption>TABLE 6. Year at which compatible fossil fuel emissions reduce to zero, where applicable, calculated from the fossil fuel emission rate smoothed with a Savitzky–Golay filter with an 11-yr window.</caption><thead><tr><th rowspan="2"></th><th colspan="8">SSP</th></tr><tr><th>1-1.9</th><th>1-2.6</th><th>2-4.5</th><th>3-7.0</th><th>4-3.4</th><th>4-6.0</th><th>5-3.4over</th><th>5-8.5</th></tr></thead><tbody><tr><td>IAM</td><td>2056</td><td>2076</td><td>—</td><td>—</td><td>2084</td><td>—</td><td>2068</td><td>—</td></tr><tr><td>ACCESS-ESM1.5</td><td>—</td><td>2076</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>CanESM5</td><td>2063</td><td>2080</td><td>—</td><td>—</td><td>2083</td><td>—</td><td>2068</td><td>—</td></tr><tr><td>CESM2</td><td>—</td><td>2082</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>CNRM-ESM2.1</td><td>2071</td><td>2083</td><td>—</td><td>—</td><td>2091</td><td>—</td><td>2078</td><td>—</td></tr><tr><td>IPSL-CM6A-LR</td><td>2070</td><td>2082</td><td>—</td><td>—</td><td>2087</td><td>—</td><td>2073</td><td>—</td></tr><tr><td>MIROC-ES2L</td><td>2060</td><td>2076</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>MPI-ESM1.2-LR</td><td>—</td><td>2084</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>NorESM2-LM</td><td>—</td><td>2086</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>UKESM-0-LL</td><td>2056</td><td>2077</td><td>—</td><td>—</td><td>2080</td><td>—</td><td>2068</td><td>—</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>FIG. 11. Model mean change in global mean temperature for all SSP experiments relative to the 1850–99 mean of the historical mean. Dashed horizontal lines indicate 1.5° and 2.0°C change.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Multimodel mean net biosphere production (GtC yr$^{-1}$) under the eight SSPs.</figcaption></figure>

cases of CESM2, CNRM, IPSL, MPI, and NorESM2, land carbon storage is also impacted by fire.

The multimodel mean land carbon uptake rates (Fig. 12) increase broadly in line with increasing atmospheric CO₂ (Fig. 1a), although the differing extents to which natural land is replaced by crop and pasture across the scenarios (Figs. 1c,d) counteract this in some cases. The high CO₂ scenario SSP5-8.5 has the greatest NBP. The SSP3-7.0, however, with second highest CO₂ has the greatest loss of forest that is replaced by the greatest pasture area and the second largest crop area of all SSPs. Additionally, SSP3-7.0 has larger contributions from non-CO₂ greenhouse gases, which do not stimulate land carbon uptake, but reduce terrestrial carbon via positive climate feedback. These factors contribute to the lowering of model mean NBP under SSP3-7.0 to below that of SSP4-6.0 (although only three ESMs performed the latter, compared to all nine for SSP3-7.0). The lowest NBP is simulated in response to SSP4-3.4; it has the third lowest CO₂ pathway, but also features the highest crop area of all scenarios, as well as fourth highest pasture area, which brings its FF emissions down approximately in line with SSP1-2.6 [section 3b(2)].

Although the model mean exhibits some interscenario spread in the time series of NBP, there is very little evidence of this in some models (ACCESS, IPSL, UKESM), while NBP in other models is very scenario dependent (CanESM5, MIROC) (Fig. S4). Figures 13b and 13c resolve the change in land carbon store (Fig. 13a) into change in vegetation carbon, and soil plus litter carbon, respectively. Figures S6–S9 show the individual ESMs’ changes in land, vegetation, soil, and litter carbon, respectively. The change in vegetation carbon is responsible, primarily, for the evolution of land carbon throughout the twenty-first century, with model mean change in soil carbon differing by less than 100 GtC across scenarios, as is true for the ESMs individually with the exception of CanESM5 (Fig. S7).

The ESMs’ natural carbon cycle behavior can be characterized by the feedback parameters introduced by the original Coupled Climate–Carbon Cycle Model Intercomparison Project (C4MIP) of the mid-2000s (Friedlingstein et al. 2006). The

C4MIP project is ongoing and forms one of the model inter-comparison projects of CMIP6; it includes a set of idealized experiments designed to allow quantification of carbon-climate and carbon-concentration feedbacks. The idealized *1pctCO2* experiments, in which CO₂ rises at 1% annually to 4 times the preindustrial level, have been performed by fully coupled, bio-geochemically coupled, and radiatively coupled variants of 11 CMIP6 ESMs, including all 9 in the present study. This combination of experiments allows the carbon-concentration feedback (denoted β₉ and βO for land and ocean, respectively, in units of GtC ppm⁻¹) and the carbon-climate feedback (γ₉ and γO, GtC °C⁻¹) to be quantified. The feedback parameters help to explain the differences in the models’ responses to the SSPs forcings. The greater the β₉, the greater the land uptake of carbon per unit increase in atmospheric CO₂, representing a negative feedback of the carbon cycle-climate system, slowing the accumulation of CO₂ in the atmosphere. A larger, negative γ₉ results in a greater loss of carbon from land per unit warming (a positive feedback, increasing atmospheric CO₂). The carbon-climate and carbon-concentration feedbacks of the nine ESMs included in the present study are included in Table A1 of Arora et al. (2020) at 2 times and 4 times preindustrial CO₂.

Among all models, CanESM5 yields the highest land carbon uptake under SSP5-8.5, leading to higher diagnosed FF emissions than all other ESMs for that scenario. The reason for this is an increase in the strength of CO₂ fertilization effect in CanESM5 compared to its predecessor CanESM2; it has the third highest β₉ (1.28 GtC ppm⁻¹ at 4 times CO₂) of the 11 models in Arora et al. (2020), for which the mean is 0.97 GtC ppm⁻¹. CanESM5 does not contain an explicit representation of nutrient constraints on land photosynthesis; the strength of the CO₂ fertilization effect over land is determined by a photosynthesis downregulation parameterization. Arora and Scinocca (2016) attempt to constrain the model parameter that determines photosynthesis downregulation, using the historical carbon budget, the amplitude of the annual cycle of globally averaged CO₂ and its trend over the historical period. Yet despite this exercise, the response of land carbon uptake in CanESM5 to increasing CO₂ is larger than other models at high CO₂. Further, of
<figure><img src="image_14.png" /><figcaption>FIG. 13. Net change (GtC) in (a) land carbon store, (b) vegetation carbon, and (c) soil + litter carbon for the SSPs. Thick lines are the multimodel means; individual models are indicated by the thin lines. Calculated from the cumulative *nbp*, and change in *cVeg*, and change in *cSoil* + *cLitter* diagnostics, respectively.</figcaption></figure>

the 11 ESMs in Arora et al. (2020), CanESM5 is unusual in having a positive carbon-climate feedback parameter, representing a negative feedback, promoting further land carbon uptake in response to warming as well as in response to rising CO₂.

CNRM has a carbon-concentration feedback comparable to CanESM5 (1.36 GtC ppm⁻¹ at 4 times CO₂), contributing to its strong land carbon uptake. It differs from CanESM5 in having a fairly strong carbon-climate feedback (−83.11 GtC K⁻¹ compared with CanESM5's +15.95 GtC °C⁻¹ and the multimodel mean of −45.07 GtC °C⁻¹ at 4 times CO₂). As the planet gets progressively warmer under the high CO₂ scenarios, CNRM's positive carbon-climate feedback dampens carbon uptake in a way that does not happen in CanESM5. MIROC also features strong gains in land carbon under SSP5-8.5 and SSP3-7.0, 406, and 368 GtC, respectively, consistent with it having the third strongest carbon-concentration feedbacks (1.12 GtC ppm⁻¹ compared with the ESM mean of 0.97). It also has a stronger than average climate-carbon feedback, of −69.57 GtC °C⁻¹, compared with the ESM mean of −45.07, but MIROC has a lower than average climate sensitivity (Arora et al. 2020) leading to warming in SSP5-8.5 of only around 4.5°C. Of the nine ESMs in this study, MIROC, CNRM, and CanESM5 are the three with the largest fraction of emissions taken up by land in Fig. 4b of Arora et al. (2020).

UKESM and IPSL exhibit comparatively weak land carbon uptake, with little dependence on CO₂ pathway due to their relatively weak carbon-climate feedbacks (−38.4 and −8.67

GtC °C⁻¹, respectively) and carbon-concentration feedbacks of 0.7 and 0.62 GtC ppm⁻¹.

Since both NorESM2 and CESM2 employ the land carbon cycle model CLM5, their responses in the scenarios are very similar (Fig. S6). They see, respectively, a rise in land carbon of 323 and 305 GtC for SSP5-8.5, and, at the lower end, 189 and 203 GtC for SSP1-2.6. This is consistent with their below average land carbon-concentration feedback. These models have different atmospheric components, and CESM2 sees greater warming than NorESM2, but with much lower than average climate-carbon feedbacks, there is little additional loss of carbon in CESM2 compared to NorESM2, even for SSP5-8.5 (Figs. S6-S9).

Arora et al. (2020) discuss some of the differences between models that affect their feedback parameters. The CO₂ fertilization effect, defined as the change in GPP per unit change in atmospheric CO₂ concentration (GtC yr⁻¹ ppm⁻¹) contributes to the carbon-concentration feedback; this is higher in the three ESMs in Arora et al. (2020), which include dynamic vegetation, including MPI, and UKESM. However, in both of these models, the turnover time of soil and vegetation carbon are lower than average. CNRM on the other hand, features strong land carbon uptake due to its relatively high turnover times, rather than a particularly strong CO₂ fertilization effect. Another differentiator is the fraction of the change in GPP (relative to the preindustrial) that gets converted to net primary production (NPP); this is a form of carbon use efficiency
<figure><img src="image_15.png" /><figcaption>FIG. 14. (a) Air-to-sea flux of CO₂ for the eight SSPs. Multimodel mean for each SSP is shown as the solid line, with model spread indicated by the shaded region around it. (b) Change in ocean carbon store for the SSPs. Thick lines are model means; thin lines are individual models.</figcaption></figure>

denoted ΔCUE to highlight that it is the change in GPP above the preindustrial that is under consideration. CanESM5 has a higher ΔCUE than all other models, which contributes to its relatively high land carbon uptake. ACCESS sees the lowest gain in land carbon throughout the scenarios because it has the lowest ΔCUE as well as a weak CO₂ fertilization effect. These combine to give it the lowest carbon-concentration feedback of the nine ESMs in this study.

The idealized experiments used to determine the carbon feedback parameters featured fixed, preindustrial nitrogen deposition and fertilizer. The SSP experiments use time-varying nitrogen deposition, providing an additional driver of land carbon uptake in those ESMs that feature explicit carbon-nitrogen interactions, limiting carbon uptake when plant-available nitrogen is scarce; this is parameterized in others (Tables 1–3 and supplemental section 1). Some models have the ability to diagnose the reduction in land carbon uptake due to the scarcity of nitrogen and other nutrients within a simulation; for example, in UKESM, NPP is calculated before and after N limitation is taken into account, each time step. In the FUN model used in CLM5 within CESM2 and NorESM2, the “carbon cost” of nitrogen uptake in units of assimilable carbon is calculated; the higher the carbon cost of nitrogen, the greater the nitrogen limitation (Fisher et al. 2019). However, assessing the full impact of nitrogen limitation on land carbon stores would require spinning up the carbon-nitrogen and carbon-only ecosystems in separate simulations, followed by parallel piControl, historical, and scenario experiments, which is prohibitively expensive.

The SSPs' land-use trajectories, forcing the myriad land-use schemes within the ESMs, superimpose another layer of complexity upon the models' biophysical and biological processes, contributing to the changes in land, soil and vegetation carbon stores seen in Figs. S6–S8. The discussion on the ESMs' representation of crop harvest in section 3a(4) is instructive when considering these changes. However, as discussed in section 3a(3), having the no-Lu variants of each SSP scenario

would be invaluable in determining the exact role of land use on the land carbon sink by comparing the net change in land carbon storage between the two.

### 5) OCEAN CARBON RESPONSE

There is much more agreement between ESMs in their ocean carbon response to the scenarios than is the case for land carbon, as shown in Figs. 14a and 14b, depicting the air-to-sea flux and change in ocean carbon store, respectively; the same are shown for the models individually in Figs. S10 and S11, respectively. The air-to-sea flux averaged across a portion of the piControl was subtracted from that of the scenarios for all ESMs as in some cases the piControl ocean was not fully in equilibrium, with a small, but not negligible, flux to or from the atmosphere. The separation between scenarios is significant in all ESMs, and their cumulative uptake totals are in close agreement (Table S5). ACCESS gains slightly more carbon than all other ESMs over the historical period, which continues throughout the SSP scenarios leading to the greatest uptake at 2100, relative to 1850, of all ESMs for its four scenarios. This is consistent with it having the largest carbon-concentration feedback (0.9 GtC ppm⁻¹) of all the ESMs included in this study for which it is known [the other eight range from 0.70 to 0.78 GtC ppm⁻¹; Table A1 in Arora et al. (2020)]. However, its negative climate-carbon feedback (−23.75 GtC °C⁻¹) is also largest in magnitude of all ESMs, the others ranging from −9.38 to −22.25 GtC °C⁻¹, promoting greater carbon loss per 1°C of warming than other ESMs, offsetting partly the gain in carbon.

All models exhibit either a stabilization or a decline in ocean carbon uptake toward the end of the twenty-first century in all scenarios except for SSP3–7.0, where carbon uptake continues to increase. The reason for this is that the transient trend in ocean carbon uptake strongly depends on the trend of atmospheric CO₂ or dCO₂/dt (Tjiputra et al. 2014). The CO₂ trend in SSP3–7.0 continues to increase by 2100, whereas it stabilizes or declines in other SSPs.

The ESMs' ocean carbon response to the scenarios is therefore much less diverse than is true for land carbon, as was
the case CMIP5 (Jones et al. 2013). Although there is diversity in the nature and complexity of ecosystems included in the ocean carbon cycle models, and to an extent the same is true of the carbonate chemistry and air-to-sea flux parameterizations, they are less conflicted than the terrestrial carbon components, which need to incorporate the effects of human interaction with land.

## 6) AIRBORNE FRACTION

The cumulative airborne fraction of fossil fuel emissions for the SSPs are shown in Fig. 15 and individually in Fig. S12.

From the initial large model differences at the end of the historical simulation (Fig. 7), there is quite strong agreement among ESMs that by the second half of the twenty-first century the airborne fraction will increase for the higher emission scenarios, SSP5-8.5 and SSP3-7.0, and reduce for the lower emission SSP1-1.9 and SSP1-2.6. Indeed, for the latter two, the models are in extremely close agreement, with the model mean CAF reducing by 0.11 over the last half of the century in both scenarios. In the higher emissions scenarios, the sinks cannot keep pace with the amount of CO₂ being released, and the fraction left in the atmosphere increases. For the lower emission scenarios, the emission rate is lower than the combined sink strength, and therefore the fraction of emissions left in the atmosphere reduces.

For the intermediate scenarios, SSP4-3.4, SSP4-6.0, and SSP5-3.4over, the CAF tends to reduce, but at a lower rate than is the case for SSP1-1.9 and SSP1-2.6, with some individual models being almost constant over the last few decades. The overshoot scenario, SSP5-3.4over, sees a slightly greater decline in CAF than SSP1-1.9 and SSP1-2.6, at 0.14 over the last 50 years. The intermodel spread is such that the highest CAF is typically around 0.15 to 0.20 greater than the lowest CAF at 2100. CNRM has the lowest CAFs almost for the duration of all scenarios.

## 4. Conclusions

For nine state-of-the-art, comprehensive Earth system models, we present the fossil fuel CO₂ emission rates compatible with the CO₂ concentration pathway used to drive them over the historical period, and under the eight Shared Socioeconomic Pathway scenarios out to 2100.

The model-mean time series of simulated historical fossil fuel CO₂ emissions compares very well with data-based estimates; the mean cumulative total from 1850 to 2014 (394 ± 59 GtC) is extremely close to the estimate from the global carbon budget 2019 (400 ± 20 GtC); only two models, however, lie within the observed estimate's uncertainty range. With the exception of CNRM and ACCESS, the models that are outside of the observed uncertainty range underestimate cumulative historical emissions. A major contributing factor to this is the stabilization of observed CO₂ concentration in the 1940s that led to lower diagnosed FF emissions in the ESMs compared to reality, which saw emissions continue to rise.

Good agreement between the observed rate of historical fossil emissions and those from ESMs provides confidence in their ability to simulate the net carbon uptake by land and ocean reliably over the historical period. Good agreement between fossil emissions diagnosed from the immensely

complex, process-based ESMs and the emissions calculated by the much simpler IAMs for the future scenarios provides confidence in our understanding of the coupled carbon cycle and climate system. Indeed, agreement between the two for the SSP experiments is quite strong, when looking at both time series and the cumulative totals of compatible fossil fuel CO₂ emissions. If this were not the case it would undermine credibility in the sequence of influence from the scenarios of population change and socioeconomic variables driving the IAM, the resulting CO₂ concentration pathway, and the physical and biogeochemical changes simulated by the ESM in response to it, as well as associated climate impacts. However, they are sufficiently close that we can have confidence that the CO₂ scenarios forcing the complex models are a realistic expression of the evolving worlds described by the Shared Socioeconomic Pathway scenarios. The fact that the IAM emissions are slightly lower than the ESMs' average would imply that the land and ocean uptake in those models is perhaps lower than that of the ESMs.

The compatible emissions time series presented here provide an uncertainty range as well as context to the single CO₂ emissions rate generated by the IAM for each SSP. Further, with the year by which society should aspire to reduce global emissions to zero gaining ever more attention, the ESMs provide an uncertainty range around this event.

As with CMIP5, the contrast between the land and the ocean carbon responses to future forcing is marked; both generations of models exhibit a strong lack of agreement in the response of the land carbon sink to the diverse range of CO₂ scenarios, with much greater agreement for ocean carbon. One primary reason land models diverge more than ocean models in terms of their future projections of carbon uptake is that ocean carbon uptake is determined almost entirely by well understood physical and chemical processes; over land, carbon uptake is determined exclusively by biological processes that are not as well understood. In addition, uncertainty over land is contributed to by several other processes including land-use change and land management practices (e.g., treatment of changing crop and pasture area, fertilizer application, irrigation, harvesting of biomass), as well as differences in nutrient cycling (nitrogen, phosphorus), all of which affect the land carbon sink.

ESMs are ever-evolving, with some already including the land carbon impact of fire-vegetation interactions. The release of CO₂ and methane from thawing permafrost (as already occurs in CESM2) as well as methane from wetlands will likely be included in future ESMs, and the development of interactive methane emissions with a 3D atmospheric CH₄ tracer, in parallel with the CO₂ equivalent, offers further opportunities to investigate feedbacks between the carbon cycle, the climate and the wider Earth system. A continuing goal in the coming decade should be to better constrain the land carbon response, despite these myriad differences.

The “E2” paired historical and hist-noLu experiments provide the best estimate of land-use and land-cover change emissions, E_LUC, in the historical experiment. The seven ESMs that have provided output from both experiments to the CMIP6 archive all underestimate land-use emissions compared to the
Airborne Fraction of Fossil Fuel emissions

<figure><img src="image_16.png" /><figcaption>FIG. 15. Cumulative airborne fraction of emitted fossil fuel CO₂ from 2015 to 2100 for the eight SSPs simulated by the ESMs. Colored lines depict individual models, with the multimodel mean as the black line.</figcaption></figure>

two reference datasets from the global carbon budget 2019 used in Fig. 3, to a greater or lesser extent, despite the GCB2019 data not accounting for the loss of additional sink capacity included in the ESMs. The ESMs at the edges of the model spread in $E_{FF}$

and $E_{LUC}$ (Figs. 3b,d) suggest an inverse relationship between the two; greater land-use emissions are paired with fewer fossil fuel emissions. This relationship suggests that in the absence of the *-noLu* variant of the SSP scenarios, the magnitude of an
ESM's historical LULCC emissions from paired experiments could potentially be used, with caution, to provide a way of qualitatively interpreting its FF emissions. If, for example, an ESM yields very high (low) historical paired-experiment LULCC emissions, its fossil fuel emissions diagnosed from future scenarios could be considered excessively low (high).

Using the E2 experimental design to allow comparison of FF, LULCC and FF + LULCC emissions between models for future scenarios, as well as for the historical period, would allow us to examine how future land use, accounting for loss of additional sink capacity and regrowth on abandoned agricultural land, impact FF emissions. It would also allow comparison between models of the response of the natural carbon cycle to realistic, rather than idealized scenarios, due primarily to changing environmental conditions rather than to the direct effects of land use on land carbon, as discussed in section 3a(3) in relation to Fig. 3f.

Aerosols, well-mixed greenhouse gases and ozone all influence the evolving climate and, therefore, land and ocean carbon uptake, and fossil emissions. The role of non-CO₂ forcings on the carbon cycle, future carbon budgets, and on the spread of the land carbon response among ESMs in particular, is an area that would benefit from further attention in the coming years. As models incorporate new processes such as interactive diffuse:direct fraction of photosynthetically active radiation, online calculation of atmospheric methane concentration from emissions, interactive nitrogen deposition as a product of atmospheric chemistry processes, as well as downgrading of carbon uptake by plants due to ozone damage, this will likely become increasingly important.

ESMs can be driven with prescribed atmospheric CO₂ concentration (C driven), as in the experiments described here, or provided with time varying fossil fuel CO₂ emissions (E driven) from which the model calculates the atmospheric CO₂ concentration. In C-driven mode, the diagnosed fossil fuel flux is the residual of the top-level carbon fluxes between the atmosphere and the land and the ocean, so the uncertainty associated with every physical, biological, and biogeochemical process that contributes to the partitioning of carbon between land and ocean accumulates in the residual FF emissions. Conversely, in E-driven mode with prescribed fossil fuel CO₂ emissions, the same uncertainties that affect the partitioning between land and ocean feed into the changing atmospheric CO₂ burden. The diagnosed fossil fuel emissions, or the changing atmospheric CO₂ burden, therefore, tend to act as integrators of model spread when comparing multiple ESMs. A useful complement to the present paper would be a similar study examining the emissions-driven experiments of the CMIP6 Earth system models.

**Acknowledgments.** SKL, CDJ, and AJW were supported by the Joint U.K. BEIS/Defra Met Office Hadley Centre Climate Programme (GA01101). This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement 641816 (CRESCENDO). DML is supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement 1852977 and by the RUBISCO Scientific Focus Area (SFA), which is sponsored by

the Regional and Global Climate Modeling (RGCM) Program in the Climate and Environmental Sciences Division (CESD) of the Office of Biological and Environmental Research in the U.S. Department of Energy Office of Science. TH was supported by the Integrated Research Program for Advancing Climate Models (Grant JPMXD0717935715), through the Ministry of Education, Culture, Sports, Science, and Technology of Japan. JS and JT acknowledge support from the Research Council of Norway through projects INES (270061), KeyClim (295046), IMPOSE (294930), and COLUMBIA (275268), as well as support from the Bjerknes Centre for Climate Research through the project LOES. RS thanks H2020 CONSTRAIN under the Grant Agreement 820829 and the support of the team in charge of the CNRM-CM climate model. Supercomputing time was provided by the Météo-France/DSI supercomputing center. RS thanks Christophe Cassou for his fruitful discussion on the paper. We thank two reviewers for their time and their constructive and helpful comments.

## REFERENCES

Arora, V. K., and J. F. Scinocca, 2016: Constraining the strength of the terrestrial CO₂ fertilization effect in the Canadian Earth System Model version 4.2 (CanESM4.2). *Geosci. Model Dev.*, **9**, 2357–2376, https://doi.org/10.5194/gmd-9-2357-2016.

—, and Coauthors, 2020: Carbon-concentration and carbon-climate feedbacks in CMIP6 models, and their comparison to CMIP5 models. *Biogeosciences*, **17**, 4173–4222, https://doi.org/10.5194/bg-17-4173-2020.

Ballantyne, A. P., and Coauthors, 2015: Audit of the global carbon budget: Estimate errors and their impact on uptake uncertainty. *Biogeosciences*, **12**, 2565–2584, https://doi.org/10.5194/bg-12-2565-2015.

Bastos, A., and Coauthors, 2016: Re-evaluating the 1940s CO₂ plateau. *Biogeosciences*, **13**, 4877–4897, https://doi.org/10.5194/bg-13-4877-2016.

Bennedsen, M., E. Hillebrand, and S. J. Koopman, 2019: Trend analysis of the airborne fraction and sink rate of anthropogenically released CO₂. *Biogeosciences*, **16**, 3651–3663, https://doi.org/10.5194/bg-16-3651-2019.

Canadell, J. G., and Coauthors, 2007: Contributions to accelerating atmospheric CO₂ growth from economic activity, carbon intensity, and efficiency of natural sinks. *Proc. Natl. Acad. Sci. USA*, **104**, 18 866–18 870, https://doi.org/10.1073/pnas.0702737104.

Ciais, P., and Coauthors, 2013: Carbon and other biogeochemical cycles. *Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 465–570.

Collins, M., and Coauthors, 2013: Long-term climate change: Projections, commitments and irreversibility. *Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 1029–1136.

Collins, W. J., and Coauthors, 2017: Aerchemmip: Quantifying the effects of chemistry and aerosols in CMIP6. *Geosci. Model Dev.*, **10**, 585–607, https://doi.org/10.5194/gmd-10-585-2017.

Dellink, R., J. Chateau, E. Lanzi, and B. Magné, 2017: Long-term economic growth projections in the shared socioeconomic pathways. *Global Environ. Change*, **42**, 200–214, https://doi.org/10.1016/j.gloenvcha.2015.06.004.

Denman, K., and Coauthors, 2007: Couplings between changes in the climate system and biogeochemistry. *Climate Change*
2007: *The Physical Science Basis*, S. Solomon et al., Eds., Cambridge University Press, 499–587.

Etheridge, D. M., L. P. Steele, R. L. Langenfelds, R. J. Francey, J.-M. Barnola, and V. I. Morgan, 1996: Natural and anthropogenic changes in atmospheric CO₂ over the last 1000 years from air in Antarctic ice and firn. *J. Geophys. Res.*, **101**, 4115–4128, https://doi.org/10.1029/95JD03410.

Eyring, V., S. Bony, G. A. Meehl, C. A. Senior, B. Stevens, R. J. Stouffer, and K. E. Taylor, 2016: Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.*, **9**, 1937–1958, https://doi.org/10.5194/gmd-9-1937-2016.

Fisher, R. A., and Coauthors, 2019: Parametric controls on vegetation responses to biogeochemical forcing in the CLM5. *J. Adv. Model. Earth Syst.*, **11**, 2879–2895, https://doi.org/10.1029/2019MS001609.

Friedlingstein, P., and Coauthors, 2006: Climate-carbon cycle feedback analysis: Results from the C4MIP model intercomparison. *J. Climate*, **19**, 3337–3353, https://doi.org/10.1175/JCLI3800.1.

—, and Coauthors, 2019: Global carbon budget 2019. *Earth Syst. Sci. Data*, **11**, 1783–1838, https://doi.org/10.5194/essd-11-1783-2019.

Gasser, T., G. P. Peters, J. S. Fuglestvedt, W. J. Collins, D. T. Shindell, and P. Ciais, 2017: Accounting for the climate-carbon feedback in emission metrics. *Earth Syst. Dyn.*, **8**, 235–253, https://doi.org/10.5194/esd-8-235-2017.

—, L. Crepin, Y. Quilcaille, R. A. Houghton, P. Ciais, and M. Obersteiner, 2020: Historical CO₂ emissions from land use and land cover change and their uncertainty. *Biogeosciences*, **17**, 4075–4101, https://doi.org/10.5194/bg-17-4075-2020.

Gidden, M. J., and Coauthors, 2019: Global emissions pathways under different socioeconomic scenarios for use in CMIP6: A dataset of harmonized emissions trajectories through the end of the century. *Geosci. Model Dev.*, **12**, 1443–1475, https://doi.org/10.5194/gmd-12-1443-2019.

Gilfillan, D., G. Marland, T. Boden, and R. Andres, 2019: Global, regional, and national fossil-fuel CO₂ emissions. CDIAC at AppState, accessed 12 May 2020, https://energy.appstate.edu/research/work-areas/cdiac-appstate.

Hajima, T., and Coauthors, 2020: Development of the MIROC-ES2L Earth System Model and the evaluation of biogeochemical processes and feedbacks. *Geosci. Model Dev.*, **13**, 2197–2244, https://doi.org/10.5194/gmd-13-2197-2020.

Hansis, E., S. J. Davis, and J. Pongratz, 2015: Relevance of methodological choices for accounting of land use change carbon fluxes. *Global Biogeochem. Cycles*, **29**, 1230–1246, https://doi.org/10.1002/2014GB004997.

Haustein, K., M. R. Allen, P. M. Forster, F. E. L. Otto, D. M. Mitchell, H. D. Matthews, and D. J. Frame, 2017: A real-time global warming index. *Sci. Rep.*, **7**, 15417, https://doi.org/10.1038/s41598-017-14828-5.

Houghton, R. A., and A. A. Nassikas, 2017: Global and regional fluxes of carbon from land use and land cover change 1850–2015. *Global Biogeochem. Cycles*, **31**, 456–472, https://doi.org/10.1002/2016GB005546.

—, J. I. House, J. Pongratz, G. R. van der Werf, R. S. DeFries, M. C. Hansen, C. Le Quéré, and N. Ramankutty, 2012: Carbon emissions from land use and land-cover change. *Biogeosciences*, **9**, 5125–5142, https://doi.org/10.5194/bg-9-5125-2012.

Hurtt, G. C., and Coauthors, 2020: Harmonization of global land use change and management for the period 850–2100 (LUH2) for CMIP6. *Geosci. Model Dev.*, **13**, 5425–5464, https://doi.org/10.5194/gmd-13-5425-2020.

IPCC, 2014: *Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects*. C. B. Field et al., Eds., Cambridge University Press, 1132 pp., http://www.ipcc.ch/pdf/assessment-report/ar5/wg2/WGIIAR5-PartA_FINAL.pdf.

Jiang, L., and B. C. O'Neill, 2017: Global urbanization projections for the shared socioeconomic pathways. *Global Environ. Change*, **42**, 193–199, https://doi.org/10.1016/j.gloenvcha.2015.03.008.

Jones, C., and Coauthors, 2013: Twenty-first-century compatible CO₂ emissions and airborne fraction simulated by CMIP5 Earth system models under four representative concentration pathways. *J. Climate*, **26**, 4398–4413, https://doi.org/10.1175/JCLI-D-12-00554.1.

—, and Coauthors, 2016: C4MIP—The Coupled Climate-Carbon Cycle Model Intercomparison Project: Experiments protocol for CMIP6. *Geosci. Model Dev.*, **9**, 2853–2880, https://doi.org/10.5194/gmd-9-2853-2016.

Joos, F., R. Meyer, M. Bruno, and M. Leuenberger, 1999: The variability in the carbon sinks as reconstructed for the last 1000 years. *Geophys. Res. Lett.*, **26**, 1437–1440, https://doi.org/10.1029/1999GL000250.

Keeling, C., T. Whorf, and M. Wahlen, 1995: Interannual extremes in the rate of rise of atmospheric carbon dioxide since 1980. *Nature*, **375**, 666–670, https://doi.org/10.1038/375666a0.

Knorr, W., 2009: Is the airborne fraction of anthropogenic CO₂ emissions increasing? *Geophys. Res. Lett.*, **36**, L21710, https://doi.org/10.1029/2009GL040613.

Lauvset, S. K., J. Tjiputra, and H. Muri, 2017: Climate engineering and the ocean: Effects on biogeochemistry and primary production. *Biogeosciences*, **14**, 5675–5691, https://doi.org/10.5194/bg-14-5675-2017.

Lawrence, D. M., and Coauthors, 2016: The Land Use Model Intercomparison Project (LUMIP) contribution to CMIP6: Rationale and experimental design. *Geosci. Model Dev.*, **9**, 2973–2998, https://doi.org/10.5194/gmd-9-2973-2016.

Le Quéré, C., M. Raupach, and J. Canadell, 2009: Trends in the sources and sinks of carbon dioxide. *Nat. Geosci.*, **2**, 831–836, https://doi.org/10.1038/ngeo689.

—, and Coauthors, 2018: Global Carbon Budget 2018. *Earth Syst. Sci. Data*, **10**, 2141–2194, https://doi.org/10.5194/essd-10-2141-2018.

Lombardozzi, D. L., Y. Lu, P. J. Lawrence, D. M. Lawrence, S. Swenson, K. W. Oleson, W. R. Wieder, and E. A. Ainsworth, 2020: Simulating agriculture in the Community Land Model version 5. *J. Geophys. Res. Biogeochem.*, **125**, e2019JG005529, https://doi.org/10.1029/2019JG005529.

MacFarling Meure, C., D. Etheridge, C. Trudinger, P. Steele, R. Langenfelds, T. van Ommen, A. Smith, and J. Elkins, 2006: Law Dome CO₂, CH₄ and N₂O ice core records extended to 2000 years BP. *Geophys. Res. Lett.*, **33**, L14810, https://doi.org/10.1029/2006GL026152.

Meinshausen, M., S. C. B. Raper, and T. M. L. Wigley, 2011: Emulating coupled atmosphere-ocean and carbon cycle models with a simpler model, MAGICC6—Part 1: Model description and calibration. *Atmos. Chem. Phys.*, **11**, 1417–1456, https://doi.org/10.5194/acp-11-1417-2011.

—, and Coauthors, 2020: The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. *Geosci. Model Dev.*, **13**, 3571–3605, https://doi.org/10.5194/gmd-13-3571-2020.

Mercado, L., N. Bellouin, S. Sitch, O. Boucher, C. Huntingford, M. Wild, and P. Cox, 2009: Impact of changes in diffuse radiation on the global land carbon sink. *Nature*, **458**, 1014–1017, https://doi.org/10.1038/nature07949.
O'Neill, B. C., and Coauthors, 2016: The scenario model inter-
comparison project (ScenarioMIP) for CMIP6. Geosci. Model
Dev., 9, 3461-3482, https://doi.org/10.5194/gmd-9-3461-2016.

Pongratz, J., C. H. Reick, R. A. Houghton, and J. I. House, 2014: Terminology as a key uncertainty in net land use and land cover change carbon flux estimates. *Earth Syst. Dyn.*, **5**, 177–195, https://doi.org/10.5194/esd-5-177-2014.

Raupach, M. R., J. G. Canadell, and C. Le Quéré, 2008: Anthropogenic and biophysical contributions to increasing atmospheric CO₂ growth rate and airborne fraction. *Biogeosciences*, **5**, 1601–1613, https://doi.org/10.5194/bg-5-1601-2008.

—, and Coauthors, 2014: The declining uptake rate of atmo-
spheric CO₂ by land and ocean sinks. *Biogeosciences*, **11**,
3453–3475, https://doi.org/10.5194/bg-11-3453-2014.

Rayner, P. J., M. R. Raupach, M. Paget, P. Peylin, and E. Koffi,
2010: A new global gridded data set of CO₂ emissions from
fossil fuel combustion: Methodology and evaluation. *J. Geophys.
Res.*, **115**, D19306, https://doi.org/10.1029/2009JD013439.

Riahi, K., and Coauthors, 2017: The shared socioeconomic path-
ways and their energy, land use, and greenhouse gas emissions
implications: An overview. *Global Environ. Change*, **42**, 153–
168, https://doi.org/10.1016/j.gloenvcha.2016.05.009.

Rogelj, J., and Coauthors, 2018: Scenarios towards limiting global mean temperature increase below 1.5°C. *Nat. Climate Change*, **8**, 325–332, https://doi.org/10.1038/s41558-018-0091-3.

Rubino, M., and Coauthors, 2013: A revised 1000 year atmospheric δ13C-CO₂ record from Law Dome and South Pole, Antarctica. *J. Geophys. Res. Atmos.*, **118**, 8482–8499, https://doi.org/10.1002/jgrd.50668.

Samir, K., and W. Lutz, 2017: The human core of the shared socio-
economic pathways: Population scenarios by age, sex and level of
education for all countries to 2100. *Global Environ. Change*, **42**,
181–192, https://doi.org/10.1016/j.gloenvcha.2014.06.004.

Sitch, S., P. Cox, W. J. Collins, and C. Huntingford, 2007: Indirect radiative forcing of climate change through ozone effects on the land-carbon sink. *Nature*, **448**, 791–794, https://doi.org/10.1038/nature06059.

Smith, C. J., and Coauthors, 2020: Effective radiative forcing and adjustments in CMIP6 models. *Atmos. Chem. Phys.*, **20**, 9591–9618, https://doi.org/10.5194/acp-20-9591-2020.

Stocker, T., and Coauthors, 2013: Technical summary. *Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 33–115.

Taylor, K. E., R. J. Stouffer, and G. A. Meehl, 2012: An overview of CMIP5 and the experiment design. *Bull. Amer. Meteor. Soc.*, **93**, 485–498, https://doi.org/10.1175/BAMS-D-11-00094.1.

Thornhill, G. D., and Coauthors, 2021: Effective radiative forcing from emissions of reactive gases and aerosols—A multimodel comparison. *Atmos. Chem. Phys.*, **21**, 853–874, https://doi.org/10.5194/acp-21-853-2021.

Tjiputra, J. F., and Coauthors, 2014: Long-term surface pCO₂ trends from observations and models. *Tellus*, **66B**, 23083, https://doi.org/10.3402/tellusb.v66.23083.

—, A. Grini, and H. Lee, 2016: Impact of idealized future stratospheric aerosol injection on the large-scale ocean and land carbon cycles. *J. Geophys. Res. Biogeosci.*, **121**, 2–27, https://doi.org/10.1002/2015JG003045.

UNFCCC, 2015: Adoption of the Paris Agreement. United Nations, 32 pp., https://unfccc.int/resource/docs/2015/cop21/eng/l09r01.pdf.

Unger, N., X. Yue, and K. L. Harper, 2017: Aerosol climate change effects on land ecosystem services. *Faraday Discuss.*, **200**, 121–142, https://doi.org/10.1039/C7FD00033B.

Vuuren, D., and Coauthors, 2011: The representative concentration pathways: An overview. *Climatic Change*, **109**, 5–31, https://doi.org/10.1007/s10584-011-0148-z.

Wei, X., M. Shao, W. Gale, and L. Li, 2014: Global pattern of soil carbon losses due to the conversion of forests to agricultural land. *Sci. Rep.*, **4**, 4062, https://doi.org/10.1038/srep04062.

Yue, C., P. Ciais, S. Luyssaert, W. Li, M. Mcgrath, J. Chang, and S. Peng, 2018: Representing anthropogenic gross land use change, wood harvest, and forest age dynamics in a global vegetation model ORCHIDEE-MICT v8.4.2. *Geosci. Model Dev.*, **11**, 409–428, https://doi.org/10.5194/gmd-11-409-2018.

Zhang, Y., and Coauthors, 2019: Increased global land carbon sink due to aerosol-induced cooling. *Global Biogeochem. Cycles*, **33**, 439–457, https://doi.org/10.1029/2018GB006051.