Supplement.pdf

Global Carbon Budget 2024

Supplementary Information

S.1 Methodology Fossil Fuel CO

emissions (EFOS)

S.1.1 Cement carbonation

From the moment it is created, cement begins to absorb CO

from the atmosphere, a process known as ‘cement

carbonation’. We estimate this CO

sink, from 1931 onwards, as the average of two studies in the literature (Cao

et al., 2020; Guo et al., 2021 extended by Huang et al., 2023). The Global Cement and Concrete Association

reports a much lower carbonation rate, but this is based on the highly conservative assumption of 0% mortar

(GCCA, 2021). Modelling cement carbonation requires estimation of a large number of parameters, including

the different types of cement material in different countries, the lifetime of the structures before demolition, of

cement waste after demolition, and the volumetric properties of structures, among others (Xi et al., 2016).

Lifetime is an important parameter because demolition results in the exposure of new surfaces to the

carbonation process. The main reasons for differences between the two studies appear to be the assumed

lifetimes of cement structures and the geographic resolution, but the uncertainty bounds of the two studies

overlap.

S.1.2 Emissions embodied in goods and services

CDIAC, UNFCCC, and BP national emission statistics ‘include greenhouse gas emissions and removals taking

place within national territory and offshore areas over which the country has jurisdiction’ (Rypdal et al., 2006),

and are called territorial emission inventories. Consumption

based emission inventories allocate emissions to

products that are consumed within a country, and are conceptually calculated as the territorial emissions minus

the ‘embodied’ territorial emissions to produce exported products plus the emissions in other countries to

produce imported products (Consumption = Territorial

–

Exports + Imports). Consumption

based emission

attribution results (e.g. Davis and Caldeira, 2010) provide additional information to territorial

based emissions

that can be used to understand emission drivers (Hertwich and Peters, 2009) and quantify emission transfers by

the trade of products between countries (Peters et al., 2011a). The consumption

based emissions have the same

global total, but reflect the trade

driven movement of emissions across the Earth's surface in response to human

activities. We estimate consumption

based emissions from 1990

2020 by enumerating the global supply chain

using a global model of the economic relationships between economic sectors within and between every country

(Andrew and Peters, 2013; Peters et al., 2011b). Our analysis is based on the economic and trade data from the

Global Trade and Analysis Project (GTAP; Narayanan et al., 2015), and we make detailed estimates for the

years 1997 (GTAP version 5), 2001 (GTAP6), and 2004, 2007, 2011, and 2014 (GTAP10.0a), covering 57

sectors and 141 countries and regions. The detailed results are then extended into an annual time series from

1990 to the latest year of the Gross Domestic Product (GDP) data (2020 in this budget), using GDP data by

expenditure in current exchange rate of US dollars (USD; from the UN National Accounts main Aggregates

database; UN, 2022) and time series of trade data from GTAP (based on the methodology in Peters et al.,

2011b). We estimate the sector

level CO

emissions using the GTAP data and methodology, add the flaring and

cement emissions from our fossil CO

dataset, and then scale the national totals (excluding bunker fuels) to

match the emission estimates from the carbon budget. We do not provide a separate uncertainty estimate for the

consumption

based emissions, but based on model comparisons and sensitivity analysis, they are unlikely to be

significantly different than for the territorial emission estimates (Peters et al., 2012b).

S.1.3 Uncertainty assessment for E

FOS

We estimate the uncertainty of the global fossil CO2 emissions at ±5% (scaled down from the published ±10 %

at ±2σ to the use of ±1σ bounds reported here; Andres et al., 2012). This is consistent with a more detailed

analysis of uncertainty of ±8.4% at ±2σ (Andres et al., 2014) and at the high

end of the range of ±5

10% at ±2σ

reported by (Ballantyne et al., 2015). This includes an assessment of uncertainties in the amounts of fuel

consumed, the carbon and heat contents of fuels, and the combustion efficiency. While we consider a fixed

uncertainty of ±5% for all years, the uncertainty as a percentage of emissions is growing with time because of

the larger share of global emissions from emerging economies and developing countries (Marland et al., 2009).

Generally, emissions from mature economies with good statistical processes have an uncertainty of only a few

per cent (Marland, 2008), while emissions from strongly developing economies such as China have

uncertainties of around ±10% (for ±1σ; Gregg et al., 2008; Andres et al., 2014). Uncertainties of emissions are

likely to be mainly systematic errors related to underlying biases of energy statistics and to the accounting

method used by each country.

S.1.4 Growth rate in emissions

We report the annual growth rate in emissions for adjacent years (in percent per year) by calculating the

difference between the two years and then normalising to the emissions in the first year: (

FOS

(

+1)

FOS

(

))/

FOS

(

)×100%. We apply a leap

year adjustment where relevant to ensure valid interpretations of

annual growth rates. This affects the growth rate by about 0.3% yr

1 (1/366) and causes calculated growth rates

to go up approximately 0.3% if the first year is a leap year and down 0.3% if the second year is a leap year.

The relative growth rate of

FOS

over time periods of greater than one year can be rewritten using its logarithm

equivalent as follows:

!"#

(

!"#

)

(2)

Here we calculate relative growth rates in emissions for multi

year periods (e.g. a decade) by fitting a linear

trend to

ln(E

FOS

)

in Eq. (2), reported in percent per year.

S.1.5 Emissions projection for 2023

To gain insight on emission trends for 2023, we provide an assessment of global fossil CO

emissions,

FOS

, by

combining individual assessments of emissions for China, USA, the EU, and India (the four countries/regions

with the largest emissions), and the rest of the world.

The methods are specific to each country or region, as described in detail below.

China

: We use a regression between monthly data for each fossil fuel and cement, and annual data for

consumption of fossil fuels / production of cement to project full

year growth in fossil fuel consumption and

cement production. The monthly data for each product consists of the following:

Coal: Production data from the National Bureau of Statistics (NBS), plus net imports from the China

Customs Administration (i.e., gross supply of coal, not including inventory changes), adjusted

using monthly production data for thermal electricity, crude steel, pig iron, coke and cement from

NBS.

Oil: Production data from NBS, plus net imports from the China Customs Administration (i.e., gross

supply of oil, not including inventory changes)

Natural gas: Same as for oil

Cement: Production data from NBS

For oil, we use data for production and net imports of refined oil products rather than crude oil. This choice is

made because refined products are one step closer to actual consumption, and because crude oil can be subject

to large market

driven and strategic inventory changes that are not captured by available monthly data.

Furthermore, refinery output in 2022 was atypically low through August of that year compared to the rest of the

year, which results in very high growth figures for the 2023 data compared to what one can likely expect for the

last four months of this year. The estimate has been adjusted down by 0.8 percentage points to account for this,

corresponding to how much lower the ratio of January

August and September

December refinery output was in

2022 compared to the average for 2014

2022.

For each fuel and cement, we make a Bayesian linear regression between year

year cumulative growth in

supply (production for cement) and full

year growth in consumption (production for cement) from annual

consumption data. In the regression model, the growth rate in annual consumption (production for cement) is

modelled as a regression parameter multiplied by the cumulative year

year growth rate from the monthly

data through August of each year for past years (through 2022). We use broad Gaussian distributions centered

around 1 as priors for the ratios between annual and through

August growth rates. We then use the posteriors for

the growth rates together with cumulative monthly supply/production data through August of 2023 to produce a

posterior predictive distribution for the full

year growth rate for fossil fuel consumption / cement production in

2023.

If the growth in supply/production through August were an unbiased estimate of the full

year growth in

consumption/production, the posterior distribution for the ratio between the monthly and annual growth rates

100

would be centered around 1. However, in practice the ratios are different from 1 (in most cases below 1). This is

101

a result of various biasing factors such as uneven evolution in the first and second half of each year, inventory

102

changes that are somewhat anti

correlated with production and net imports, differences in statistical coverage,

103

and other factors that are not captured in the monthly data.

104

For fossil fuels, the mean of the posterior distribution is used as the central estimate for the growth rate in 2023,

105

while the edges of a 68% credible interval (analogous to a 1

sigma confidence interval) are used for the upper

106

and lower bounds.

107

USA

: We use emissions estimated by the U.S. Energy Information Administration (EIA) in their Short

Term

108

Energy Outlook (STEO) for emissions from fossil fuels to get both YTD and a full year projection (EIA, 2023).

109

The STEO also includes a near

term forecast based on an energy forecasting model which is updated monthly

110

(we use the November 2023 edition), and takes into account expected temperatures, household expenditures by

111

fuel type, energy markets, policies, and other effects. We combine this with our estimate of emissions from

112

cement production using the monthly U.S. cement clinker production data from USGS for January

August

113

2023, assuming changes in clinker production over the first part of the year apply throughout the year.

114

India

: We use monthly emissions estimates for India updated from Andrew (2020b) through August

October

115

2023. These estimates are derived from many official monthly energy and other activity data sources to produce

116

direct estimates of national CO

emissions, without the use of proxies. Emissions from coal are then extended to

117

October using a regression relationship based on power generated from coal, coal dispatches by Coal India Ltd.,

118

the composite PMI, time, and days per month. For the last 3

5 months of the year, each series is extrapolated

119

assuming typical (pre

2019) trends.

120

: We use a refinement to the methods presented by Andrew (2021), deriving emissions from monthly energy

121

data reported by Eurostat. Some data gaps are filled using data from the Joint Organisations Data Initiative

122

(JODI, 2022). Sub

annual cement and cement

clinker production data are limited, but data for Germany, Poland

123

and Spain, the three largest producers, suggest a decline of over 8%. For fossil fuels this provides estimates

124

through July

September, varying by fuel. We extend coal emissions through October using a regression model

125

built from generation of power from hard coal, power from brown coal, and the number of working days in

126

Germany, the biggest coal consumer in the EU. These are then extended through the end of the year assuming

127

typical trends. We extend oil emissions by building a regression model between our monthly CO

estimates and

128

oil consumption reported by the EIA for Europe in its Short

Term Energy Outlook (November edition), and then

129

using this model with EIA’s monthly forecasts. For natural gas, the strong seasonal signal allows the use of the

130

bias

adjusted Holt

Winters exponential smoothing method (Chatfield, 1978), although this comes with larger

131

uncertainty given the unusual energy situation in Europe in 2022

23.

132

Rest of the world

: We use the close relationship between the growth in GDP and the growth in emissions

133

(Raupach et al., 2007) to project emissions for the current year. This is based on a simplified Kaya Identity,

134

whereby E

FOS

(GtC yr

) is decomposed by the product of GDP (USD yr

) and the fossil fuel carbon intensity of

135

the economy (I

FOS

; GtC USD

) as follows:

136

퐸

*+,

퐺퐷푃

퐼

*+,

(3)

137

Taking a time derivative of Equation (3) and rearranging gives:

138

!"#

-./

#-./

!"#

(4)

139

where the left

hand term is the relative growth rate of E

FOS

, and the right

hand terms are the relative growth

140

rates of GDP and I

FOS

, respectively, which can simply be added linearly to give the overall growth rate.

141

The I

FOS

is based on GDP in constant PPP (Purchasing Power Parity) from the International Energy Agency

142

(IEA) up to 2017 (IEA/OECD, 2019) and extended using the International Monetary Fund (IMF) growth rates

143

through 2022 (IMF, 2023). Interannual variability in I

FOS

is the largest source of uncertainty in the GDP

based

144

emissions projections. We thus use the standard deviation of the annual IFOS for the period 2013

2022 as a

145

measure of uncertainty, reflecting a ±1σ as in the rest of the carbon budget. For rest

world oil emissions

146

growth, we use the global oil demand forecast published by the EIA less our projections for the other four

147

regions, and estimate uncertainty as the maximum absolute difference over the period available for such

148

forecasts using the specific monthly edition (e.g. August) compared to the first estimate based on more solid

149

data in the following year (April).

150

Bunkers

: Given the divergence in behaviour of international shipping from countries’ emissions since the

151

COVID

19 pandemic, we project international bunkers separately using sub

annual data on international

152

aviation from the OECD (Clarke et al., 2022) and international shipping from MarineBenchmark and IMF

153

(Cerdeiro et al., 2020).

154

World

: The global total is the sum of each of the countries and regions.

155

156

S.2 Methodology CO

emissions from land

use, land

use change and forestry (E

LUC

)

157

The net CO

flux from land

use, land

use change and forestry (E

LUC

, called land

use change emissions in the

158

rest of the text) includes CO

fluxes from deforestation, afforestation, logging and forest degradation (including

159

harvest activity), shifting cultivation (cycle of cutting forest for agriculture, then abandoning), regrowth of

160

forests following wood harvest or abandonment of agriculture, peat burning, and peat drainage. Land

161

management activities are only partly included in our land

use change emissions estimates (Table S1). Some

162

land

use change and land

management activities cause emissions of CO

to the atmosphere, while others

163

remove CO

from the atmosphere. E

LUC

is the net sum of emissions and removals due to all anthropogenic

164

activities considered. Our annual estimates for 1960

2022 are provided as the average of results from four

165

bookkeeping approaches (Supplement S.2.1 below): the Bookkeeping of Land Use Emissions model (BLUE;

166

Hansis et al., 2015), the compact Earth system model OSCAR (Gasser et al., 2020), an estimate from Houghton

167

and Castanho (2023; hereafter H&C2023), and the Land

Use Change Emissions model (LUCE; Qin et al.,

168

2024). Peat emissions are added from external datasets (see Supplement S.2.1 below). BLUE and OSCAR are

169

updated with new land

use forcing data covering the time period until 2023. All four data sets are extrapolated

170

to provide a projection for 2024 (see Supplement S.2.5 below). In addition, we use results from Dynamic Global

171

Vegetation Models (DGVMs; see Supplement S.2.2 and Table 4) to help quantify the uncertainty in E

LUC

172

(Supplement S.2.4), and thus better characterise the robustness of annual estimates and trends. In this budget, we

173

follow the scientific E

LUC

definition as used by global carbon cycle models, which counts fluxes due to

174

environmental changes on managed land towards S

LAND

, as opposed to the national greenhouse gas inventories

175

(NGHGIs) under the UNFCCC, most of which include them in E

LUC

and thus often report smaller land

use

176

emissions (Grassi et al., 2018; Petrescu et al., 2020). Following the methodology of Grassi et al. (2023), we

177

provide harmonised estimates of the two approaches further below (see Supplement S.2.3).

178

S.2.1 Bookkeeping models

179

emissions and removals from land

use change are calculated by four bookkeeping models. These are based

180

on the original bookkeeping approach of Houghton (2003), which keeps track of the carbon stored in vegetation

181

and soils before and after a land

use change event (transitions between various natural vegetation types,

182

croplands, and pastures). Literature

based response curves describe the decay of vegetation and soil carbon,

183

including carbon transfer to product pools of different lifetimes, as well as carbon uptake due to regrowth. In

184

addition, the bookkeeping models represent long

term degradation of primary forest as lowered standing

185

vegetation and soil carbon stocks in secondary forests, and include forest management practices such as wood

186

harvests.

187

BLUE, LUCE and H&C2023 exclude the transient response of land ecosystems to changes in climate,

188

atmospheric CO

, and other environmental factors, and base the carbon densities of soil and vegetation on

189

contemporary data from literature and inventory data. Since carbon densities thus remain fixed over time, the

190

additional sink capacity that ecosystems provide in response to CO

fertilisation and some other environmental

191

changes are not captured by these models (Pongratz et al., 2014). OSCAR includes this transient response, and it

192

follows a theoretical framework (Gasser and Ciais, 2013) that allows separating bookkeeping land

use

193

emissions and the loss of additional sink capacity. Only the former is included here, while the latter is discussed

194

in Supplement S.6.4. The bookkeeping models differ in (1) computational units (spatially explicit treatment of

195

land

use change at 0.25° resolution for BLUE and LUCE, country

level for H&C2023 and OSCAR), (2)

196

processes represented (see Table S1), and (3) carbon densities assigned to vegetation and soils for different

197

types of vegetation (literature

based for BLUE and H&C2023, calibrated to DGVMs for OSCAR, mainly

198

literature

based but additionally considering the impact of land cohort age on secondary land carbon stocks for

199

LUCE). A notable difference between models exists with respect to the treatment of shifting cultivation:

200

H&C2023 assumes that forest loss

—

derived from the Global Forest Resources Assessment (FRA; FAO,

201

2020)

—

in excess of increases in cropland and pastures

—

derived from FAOSTAT (FAO, 2021)

—

represents an

202

increase in shifting cultivation. If the excess loss of forests in a year is negative, it is assumed that shifting

203

cultivation is returned to forest. Historical areas in shifting cultivation are defined taking into account country

204

based estimates of areas in fallow in 1980 (FAO/UNEP, 1981) and expert opinion (from Heinimann et al.,

205

2017). In contrast, BLUE, OSCAR, and LUCE include subgrid

scale transitions between all vegetation types.

206

Furthermore, H&C2023 assumes conversion of natural grasslands to pasture, while BLUE, OSCAR, and LUCE

207

allocate pasture transitions proportionally to all natural vegetation that exists in a grid

cell. This is one reason

208

for generally higher emissions in BLUE and OSCAR. In this GCB, we split CO

emissions into emissions from

209

permanent deforestation and from deforestation for shifting cultivation. Similarly, we separate the forest (re

210

)growth estimates into (re

)growth from af/reforestation and from regrowth associated with shifting cultivation.

211

This distinction is insightful with regard to the levers on the reduction of net emissions: as deforestation for

212

shifting cultivation is only temporary, the associated CO

emissions cannot easily be avoided without

213

compromising the CO

removals from regrowth in shifting cultivation cycles. By contrast, permanent

214

deforestation is typically not directly related to af/reforestation. Stopping deforestation for permanent

215

agricultural expansion and increasing the forest area provide two independent levers for net emissions reduction.

216

Bookkeeping models do not directly capture carbon emissions from the organic layers of drained peat soils nor

217

from peat fires. Particularly the latter can create large emissions and interannual variability due to synergies of

218

land

use and climate variability in equatorial Southeast Asia, especially during El

Niño events. We add peat fire

219

emissions based on the Global Fire Emission Database (GFED4s; van der Werf et al., 2017) to the bookkeeping

220

models’ output. Peat fire emissions are calculated by multiplying the mass of dry matter emitted by peat fires

221

with the C emission factor for peat fires indicated in the GFED4s database. Emissions from deforestation and

222

degradation fires used for extrapolating the H&C2023 data beyond 2020 and to derive the 2023 projection of all

223

three models (see below) are calculated analogously. The satellite

derived GFED4s estimates of peat fire

224

emissions start in 1997. For the previous years, we follow the approach by Houghton and Nassikas (2017),

225

which linearly ramps up from zero emissions in 1980 to 0.04 GtC yr

in 1996, reflecting the onset of major

226

clearing of peatlands in equatorial Southeast Asia in the 1980s.

227

We further add estimates of peat drainage emissions, combining estimates from three spatially explicit datasets.

228

We employ FAO peat drainage emissions 1990

–

2022 from croplands and grasslands (Conchedda and Tubiello,

229

2020; FAO, 2023), peat drainage emissions 1700

–

2010 from simulations with the DGVM ORCHIDEE

PEAT

230

(Qiu et al., 2021), and peat drainage emissions 1701

–

2023 from simulations with the DGVM LPX

Bern v1.5

231

(Lienert and Joos, 2018; Müller and Joos, 2021), the latter applying the updated LUH2

GCB2024 forcing as

232

also used by BLUE, OSCAR, LUCE, and the DGVMs. The LPX

Bern simulations started from a transient run

233

over the last deglaciation (

20,050 to 1700 AD) following Müller and Joos (2020) and are forced by changes in

234

climate, atmospheric CO

, nitrogen deposition/input, and land

use changes. Simulations were done with and

235

without prescribing land

use changes since 1700 AD. The difference between the simulations represents

236

anthropogenic peat drainage emissions. To account for internal variability, we used the median peat drainage

237

emissions from a 20

member ensemble. In LPX

Bern, peat carbon is stored in (i) active peatlands, (ii) former

238

peatlands (“natural”), and (iii) former peatlands under anthropogenic use. We average the two CO

emission

239

cases from Müller and Joos (2021), assuming that half the peat carbon is lost immediately to the atmosphere

240

after transformation from active to former peatland, while the rest decays slowly, pending on local temperature

241

and soil moisture. The LPX

Bern peat drainage emissions show a very high emission peak in Russia in 1959

242

followed by very low emissions in 1960. This peak can be attributed to an artefact in the HYDE3.4 dataset,

243

which was corrected for Brazil and the Democratic Republic of the Congo in GCB2022 (Friedlingstein et al.

244

2022b) but remains for Russia where it strongly impacts the LPX

Bern peat drainage estimates in 1959 and

245

1960. To correct for this unrealistic peak, we replace the LPX

Bern peat drainage emissions in Russia in 1959

246

and 1960 by the average of the estimates in 1958 and 1961. FAO data are extrapolated to 1850

2023 by keeping

247

the post

2020 emissions constant at 2020 levels and by linearly increasing tropical peat drainage emissions

248

between 1980 and 1990 starting from 0 GtC yr

in 1980 (consistent with H&N2017’s assumption, Houghton

249

and Nassikas, 2017), and by keeping pre

1990 emissions from the often old drained areas of the extra

tropics

250

constant at 1990 emission levels. ORCHIDEE

PEAT data are extrapolated to 2011

2023 by replicating the

251

average emissions in 2000

2010 (pers. comm. C. Qiu). Further, ORCHIDEE

PEAT only provides peat drainage

252

emissions north of 30°N, and thus we fill the regions south of 30°N by the average peat drainage emissions from

253

FAO and LPX

Bern. The final peat drainage emissions are calculated as the average of the estimates from the

254

three different peat drainage datasets. The net E

LUC

values indicated in the manuscript are the sum of E

LUC

255

estimates from bookkeeping models, peat fire emissions, and peat drainage emissions.

256

The four bookkeeping estimates used in this study differ with respect to the land

use change data used to drive

257

the models. H&C2023 base their estimates directly on the Forest Resource Assessment (FRA) of FAO, which

258

provides statistics on forest

area change and management at intervals of five years currently updated until 2020

259

(FAO, 2020). The data is based on country reporting to FAO and may include remote

sensing information in

260

more recent assessments. Changes in land use other than forests are based on annual, national changes in

261

cropland and pasture areas reported by FAO (FAO, 2021). BLUE and LUCE use the harmonised land

use

262

change data LUH2

GCB2024 covering the period 850

2023 (an update to the previously released LUH2 v2h

263

dataset; Hurtt et al., 2017; Hurtt et al., 2020), which was also used as input to the DGVMs (Supplement S.2.2).

264

LUH2

GCB2024 provides land

use change data at 0.25° spatial resolution based on the FAO data (as described

265

in Supplement S.2.2) as well as the HYDE3.4 dataset (Klein Goldewijk et al., 2017a, 2017b), considering

266

subgrid

scale transitions between primary forest, secondary forest, primary non

forest, secondary non

forest,

267

cropland, pasture, rangeland, and urban land (Hurtt et al., 2020; Chini et al., 2021). LUH2

GCB2024 provides a

268

distinction between rangelands and pasture, based on inputs from HYDE. Rangeland establishment in forests is

269

assumed to transform forests to grasslands, rangeland establishment in non

forest primary vegetation degrades

270

primary to secondary vegetation, and rangeland establishment in non

forest secondary vegetation has no effect

271

(e.g., browsing on shrubland) (Ma et al., 2020). This case distinction is implemented in BLUE based on a forest

272

mask provided with LUH2

GCB2021. OSCAR was run with both LUH2

GCB2024 and FAO/FRA, where the

273

drivers of the latter were linearly extrapolated to 2023 using their 2015

2020 trends. The best

guess OSCAR

274

estimate used in our study is a combination of results for LUH2

GCB2024 and FAO/FRA land

use data and a

275

large number of perturbed parameter simulations weighted against a constraint (the cumulative S

LAND

over

276

1960

2022 of last year’s GCB). As the record of H&C2023 ends in 2020, we extend it up to 2023 by adding the

277

yearly anomalies of the emissions from tropical deforestation and degradation fires from GFED4s between 2020

278

and 2022 to the model’s estimate for 2020 (emissions from peat fires and peat drainage are added to all models

279

later in the process).

280

The annual E

LUC

from 1850 onwards is calculated as the average of the estimates from BLUE, H&C2023,

281

OSCAR, and LUCE. For the cumulative numbers starting in 1750, emission estimates between 1750

1850 are

282

added based on the average of four earlier publications (30 ± 20 GtC 1750

1850, rounded to nearest 5; Le Quéré

283

et al., 2016).

284

We provide a split of net E

LUC

into component fluxes to better identify reasons for divergence between

285

bookkeeping estimates and to give more insight into the drivers of net E

LUC

. This split distinguishes between

286

emissions from deforestation (including due to shifting cultivation), removals from forest (re

)growth (including

287

regrowth in shifting cultivation cycles), fluxes from wood harvest and other forest management (i.e., emissions

288

in forests from slash decay and emissions from product decay following wood harvesting, removals from

289

regrowth after wood harvesting, and fire suppression), emissions from peat drainage and peat fires, and

290

emissions and removals associated with all other land

use transitions. Additionally, we split deforestation

291

emissions into emissions from permanent deforestation and emissions from deforestation in shifting cultivation

292

cycles, and we split removals from forest (re

)growth into forest (re

)growth due to af/reforestation and forest

293

regrowth in shifting cultivation cycles. This split helps to identify the emission reductions that would be

294

achievable by halting permanent deforestation, and the removals that are caused by permanently increasing the

295

forest cover through re/afforestation. Forest (re

)growth due to af/reforestation is calculated using a slightly

296

updated method compared to GCB2023, now following the method used to calculate CDR due to

297

re/afforestation in the 2nd State of CDR Report (Pongratz et al., 2024). E

LUC

data are provided as global sums,

298

as spatially explicit estimates at 0.25° spatial resolution (i.e., the native LUH2 resolution), and for 199 countries

299

(based on the list of UNFCCC parties). Spatially explicit E

LUC

estimates for BLUE and LUCE are directly

300

available at 0.25°. For OSCAR and H&C2023, the country

level estimates were scaled to 0.25° based on the

301

patterns of gross emissions and gross removals in BLUE (see Schwingshackl et al. 2022 for more details about

302

the methodology). The gridded net E

LUC

estimates of BLUE, LUCE, OSCAR, and H&C2023 are averaged, and

303

the gridded estimates of peat drainage emissions (average of FAO, LPX

Bern, and ORCHIDEE

PEAT) and of

304

peat fire emissions (from GFED4s) are added. Country

level estimates for the gridded datasets (BLUE, LUCE,

305

LPX

Bern, ORCHIDEE

PEAT, GFED4s) are calculated based on a country map from Eurostat (Eurostat,

306

2024), which was remapped to 0.25°. In case multiple countries are present in a 0.25° grid cell, the E

LUC

307

estimates are allocated proportional to each country’s land fraction in that grid cell.

308

309

S.2.2 Dynamic Global Vegetation Models (DGVMs)

310

Land

use change CO

emissions are also estimated by an ensemble of 20 DGVMs. The DGVMs account for

311

deforestation and regrowth, the most important components of E

LUC

, but they do not represent all processes

312

resulting directly from human activities on land (Table S1). All DGVMs represent processes of vegetation

313

growth and mortality, as well as decomposition of dead organic matter associated with natural cycles, and

314

include the vegetation and soil carbon response to increasing atmospheric CO

concentration, to climate

315

variability and to climate change. Most models explicitly simulate the coupling of carbon and nitrogen cycles

316

and account for atmospheric N deposition and N fertilisers (Table S1). The DGVMs are independent from the

317

other budget terms except for their use of atmospheric CO

concentration to calculate the fertilisation effect of

318

on plant photosynthesis.

319

All DGVMs use the LUH2

GCB2024 dataset as input, which includes the HYDE cropland/grazing land dataset

320

(Klein Goldewijk et al., 2017a, 2017b), and some additional information on land

use transitions, land

use

321

management activities and wood harvest. This includes annual, quarter

degree (regridded from 5 minute

322

resolution), fractional data on cropland and pasture from HYDE3.4.

323

DGVMs that do not simulate subgrid

scale transitions (i.e., those estimating net land

use emissions; see Table

324

S1) used the HYDE information on agricultural area change. For all countries, with the exception of Brazil, the

325

Democratic Republic of the Congo, Indonesia, and China these data are based on the available annual FAO

326

statistics of change in agricultural land area available from 1961 up to and including 2017. The FAO

327

retrospectively revised their reporting for the Democratic Republic of the Congo, which was newly available

328

until 2020 as reported in GCB2022. In addition to FAO country

level statistics, the HYDE3.4 cropland/grazing

329

land dataset is constrained spatially based on multi

year satellite land cover maps from ESA CCI LC (see

330

below). The extension of HYDE beyond the years that were directly informed by data was done as part of the

331

LUH2 methodology this year and was a simple extension of the previous 5

year trend. The actual years for this

332

extension varied by country since some countries were based on FAO data (2021), some used the China data

333

(2019), and some used MapBiomas data (Brazil and Indonesia, 2022). This methodology is not appropriate for

334

countries that have experienced recent rapid changes in the rate of land

use change, e.g. Brazil which has

335

experienced a recent upturn in deforestation. For Brazil and Indonesia we replace FAO state

level data for

336

cropland and grazing land in HYDE by those from the satellite

based land cover dataset MapBiomas (collection

337

7) for 1985

2022 (Brazil) (Souza et al. 2020) and 2000

2022 (Indonesia). ESA

CCI is used to spatially

338

disaggregate as described below.. The pre

1985 period is scaled with the per capita numbers from 1985 from

339

MapBiomas, so this transition is smooth.

340

HYDE uses satellite imagery from ESA

CCI from 1992

2018 for more detailed yearly allocation of cropland

341

and grazing land, with the ESA area data scaled to match the FAO annual totals at country

level. The original

342

300 metre spatial resolution data from ESA was aggregated to a 5 arc minute resolution according to the

343

classification scheme as described in Klein Goldewijk et al. (2017a).

344

DGVMs that simulate subgrid

scale transitions (i.e., those estimating gross land

use emissions; see Table S1)

345

use more detailed land use transition and wood harvest information from the LUH2

GCB2024 data set. LUH2

346

GCB2024 is an update of the comprehensive harmonised land

use data set (Hurtt et al., 2020), that includes

347

fractional data on primary and secondary forest vegetation, as well as all underlying transitions between land

348

use states (850

2023; Hurtt et al., 2011, 2017, 2020; Chini et al., 2021; Table S1). This data set consists of 0.25°

349

fractional areas of land

use states and all transitions between those states, including a new wood harvest

350

reconstruction, new representation of shifting cultivation, crop rotations, management information including

351

irrigation and fertiliser application. The land

use states include five different crop types in addition to splitting

352

grazing land into managed pasture and rangeland. Wood harvest patterns are constrained with Landsat

based

353

tree cover loss data (Hansen et al. 2013). Updates of LUH2

GCB2024 over last year’s version (LUH2

354

GCB2023) are using the most recent HYDE release. HYDE4.3 is based on new FAO inputs for years 1961

355

2021, new MapBiomas inputs for Brazil (for years 1985

2022) and Indonesia (for years 2000

2022) and new

356

cropland data for China from Yu et al. 2022 (for years 1900

2019).

357

We use updated FAO wood harvest data for all dataset years from 1961 to 2022, and linearly extended to the

358

year 2023. The HYDE3.4 population data is also used to extend the wood harvest time series back in time.

359

Other wood harvest inputs (for years prior to 1961) remain the same in LUH2. These updates in the land

use

360

forcing are shown in Figure S7 in comparison to LUH2

GCB2022 and LUH2

GCB2023. DGVMs implement

361

land

use change in different ways (e.g. an increased cropland fraction in a grid cell can either be at the expense

362

of grassland, shrubs, or forest, the latter resulting in deforestation; land cover fractions of the non

agricultural

363

land differ between models). Similarly, model

specific assumptions are applied to convert deforested biomass or

364

deforested area, and other forest product pools into carbon, and different choices are made regarding the

365

allocation of rangelands as natural vegetation or pastures.

366

The difference between two DGVMs simulations (see Supplement S.4.1 below), one forced with historical

367

changes in land

use and a second one with time

invariant pre

industrial land cover and pre

industrial wood

368

harvest rates, allows quantification of the dynamic evolution of vegetation biomass and soil carbon pools in

369

response to land

use change in each model (E

LUC

). Using the difference between these two DGVM simulations

370

to diagnose E

LUC

means the DGVM estimate includes the loss of additional sink capacity (around 0.4 ± 0.3 GtC

371

1; see Section 2.10 and Supplement S.6.4), while the bookkeeping model estimate does not.

372

As a criterion for inclusion in this carbon budget, we only retain models that simulate a positive E

LUC

during the

373

1990s, as assessed in the IPCC AR4 (Denman et al., 2007) and AR5 (Ciais et al., 2013). All DGVMs met this

374

criterion.

375

376

S.2.3 Translation between NGHGIs and E

LUC

377

Land

use emissions estimates from bookkeeping models and from national GHG Inventories (NGHGIs) show a

378

large gap (see Figure 8 and Table S10). This gap is due to different approaches for calculating “anthropogenic”

379

fluxes related to land

use change and land management (Grassi et al. 2018). Land sinks due to

380

environmental change on managed lands are treated as non

anthropogenic in the global carbon budget, while

381

they are generally considered as anthropogenic in NGHGIs (“indirect anthropogenic fluxes”; Eggleston et al.,

382

2006). Building on previous studies (Grassi et al. 2021), we implement an approach that adds the DGVM

383

estimates of CO

fluxes due to environmental change from managed forest areas (part of S

LAND

) to the E

LUC

384

estimate from bookkeeping models. This sum is expected to be conceptually more comparable to NGHGI

385

estimates than E

LUC

386

LUC

data are taken from bookkeeping models, in line with the global carbon budget approach. To determine

387

LAND

in managed forest, the following steps were taken: Spatially gridded data of “natural” forest NBP (S

LAND

388

i.e., including carbon fluxes due to environmental change and excluding land use change fluxes) were obtained

389

from DGVMs using S2 runs from the TRENDY v13 dataset. Results were first masked with a forest map that is

390

based on tree cover data from Hansen et al. (2013). To perform the conversion “tree” cover to “forest” cover, we

391

exclude gridcells with less than 20% tree cover and isolated pixels with maximum connectivity less than 0.5 ha

392

following the FAO definition of forest. Forest NBP is then further masked with a map of “intact” forest for the

393

year 2013, i.e. forest areas characterised by no remotely detected signs of human activity (Potapov et al. 2017).

394

This way, we obtained S

LAND

in “intact” and “non

intact” forest areas, which previous studies (Grassi et al.

395

2021) indicated to be a good proxy, respectively, for “unmanaged” and “managed” forest areas in the NGHGI.

396

Note that only a subset of models had forest NBP at grid cell level. For the other DGVMs, when a grid cell had

397

forest, all the NBP in that grid cell was allocated to forest. Since S2 simulations use pre

industrial forest cover

398

masks that are at least 20% larger than today’s forest (Hurtt et al. 2020), we corrected this NBP by a ratio

399

between observed (based on Hansen et al. 2013) and prescribed (from DGVMs) forest cover. This ratio is

400

calculated for each individual DGVM that provides information on prescribed forest cover, and a common ratio

401

(median ratio of this subset of models) is used. The details of the method used are explained in a GitHub

402

repository (Alkama, 2022).

403

LULUCF data from NGHGIs are from Grassi et al. (2023), updated up to August 2024. While Annex I countries

404

report a complete time series 1990

2021, gap

filling was applied for Non

Annex I countries through linear

405

interpolation between two points and/or through extrapolation backward (till 2000) and forward (till 2021) using

406

the single closest available data. For all countries, the estimates of the years 2022 and 2023 are assumed to be

407

equal to those of 2021. The managed forest area, used to filter SLAND data from DGVMs to derive the natural

408

land sink in managed forests, accounts for temporal dynamics from 2000 to 2023. This data includes all CO

409

fluxes from land considered managed, which in principle encompasses all land uses (forest land, cropland,

410

grassland, wetlands, settlements, and other land), changes among them, emissions from organic soils (i.e., from

411

peat drainage) and from fires. In practice, although almost all Annex I countries report all land uses, many non

412

Annex I countries report only on deforestation and forest land, and only few countries report on other land uses.

413

In most cases, NGHGIs include most of the natural response to recent environmental change because they use

414

direct observations (e.g., national forest inventories) that do not allow separating direct and indirect

415

anthropogenic effects (Eggleston et al., 2006).

416

Last, we also used the gridded data of net land flux from 14 atmospheric inversion systems (Table S4) to get an

417

additional estimate of land

use fluxes in managed land. We applied a correction for riverine transport (see

418

Supplement S.5.1.) and multiplied the resulting values with the fraction of managed land in each grid cell for

419

each inversion. For this purpose, we used masks of managed land from Grassi et al. (2023) available for the

420

years 1994, 2002, 2010, and 2016. We linearly interpolated the masks in time and replicated the 2016 mask in

421

the years 2017

2023. Subsequently, we applied another correction for lateral transport due to international wood

422

and crop trade (data from

Deng et al. 2024

). The obtained values are summed globally and compared to the

423

NGHGI estimates and the translated E

LUC

estimates.

424

Figure 8 and Table S10 shows the resulting translation of global carbon cycle models' land flux definitions to

425

that of the NGHGI (discussed in Section 3.2.2). For comparison we also show LULUCF estimates from

426

FAOSTAT (FAO, 2024), which include emissions from net forest conversion and fluxes on forest land

427

(Tubiello et al., 2021) as well as CO

emissions from peat drainage and peat fires. Forest land stock change data

428

for 2021

2023 are carried forward from the 2020 estimates. The FAO data shows global emissions of 0.30 GtC

429

averaged over 2014

2023, in contrast to the removals of

0.76 GtC yr

estimated by the gap

filled NGHGI

430

data. Most of this difference is attributable to different scopes: a focus on carbon fluxes for the NGHGI and a

431

focus on land

use area and biomass estimates for FAO. In particular, the NGHGI data includes a larger forest

432

sink for non

Annex 1 countries resulting from a more complete coverage of non

biomass carbon pools and non

433

forest land uses. NGHGI and FAO data also differ in terms of underlying data on forest land (Grassi et al.,

434

2022).

435

436

S.2.4 Uncertainty assessment for E

LUC

437

Differences between the bookkeeping models and DGVMs originate from three main sources: different

438

methodologies, which among others lead to inclusion of the loss of additional sink capacity in DGVMs (see

439

Supplement S.6.4), different underlying land

use/land cover datasets, and different processes represented (Table

440

S1). We examine both the results from DGVMs and from the bookkeeping method and use the resulting

441

variations as a way to characterise the uncertainty in E

LUC

442

Despite the existing differences, the E

LUC

estimate from the DGVM multi

model mean is consistent with the

443

average of the emissions from the bookkeeping models (Table 5). However, there are large differences among

444

individual DGVMs (standard deviation at 0.6 GtC yr

; Table 5), between the bookkeeping estimates (standard

445

deviation at 0.3 GtC yr

for cumulative emissions in 1850

2022), and between the H&C2023 model and its

446

previous model version H&N2017 (average difference 1850

2015 of 0.2 GtC yr

; see Table 1 in Houghton and

447

Castanho, 2023). A factorial analysis of differences between BLUE and H&N2017 (the precursor of H&C2023)

448

attributed them particularly to differences in carbon densities between primary and secondary vegetation (Bastos

449

et al., 2021). Earlier studies additionally showed the relevance of the different land

use forcing as applied (in

450

updated versions) also in the current study (Gasser et al., 2020). Ganzenmüller et al. (2022) showed that E

LUC

451

estimates with BLUE are substantially smaller when the model is driven by a new high

resolution land

use

452

dataset (HILDA+). They identified shifting cultivation and the way it is implemented in LUH2 as a main reason

453

for this divergence. They further showed that a higher spatial resolution reduces the estimates of both gross

454

emissions and gross removals because successive transitions are not adequately represented at coarser

455

resolution, which has the effect that

—

despite capturing the same extent of transition areas

—

overall less area

456

remains pristine at the coarser compared to the higher resolution.

457

The uncertainty in E

LUC

of ±0.7 GtC yr

reflects our best value judgement that there is at least 68% chance

458

(±1σ) that the true land

use change emissions lie within the given range, for the range of processes considered

459

here. Prior to the year 1959, the uncertainty in E

LUC

is taken from the standard deviation of the DGVMs. We

460

assign low confidence to the annual estimates of E

LUC

because of the inconsistencies among estimates and

461

because of the difficulties to quantify some of the processes with DGVMs.

462

463

S.2.5 Land

use emissions projection for 2024

464

We project the 2024 land

use emissions for BLUE, H&C2023, OSCAR, and LUCE based on their E

LUC

465

estimates for 2023 and on the interannual variability of peat fires and tropical deforestation and degradation fires

466

as estimated using active fire data (MCD14ML; Giglio et al., 2016). The latter scales almost linearly with GFED

467

emissions estimates over large areas (van der Werf et al., 2017), and thus allows for tracking fire emissions in

468

deforestation and tropical peat zones in near

real time. Peat drainage is assumed to be unaltered, as it has low

469

interannual variability. We project the 2024 land

use emissions for BLUE, H&C2023, OSCAR, and LUCE

470

based on their E

LUC

estimates for 2023 and add the change in carbon emissions from peat fires and tropical

471

deforestation and degradation fires (2024 emissions relative to 2023 emissions) from GFED4s. The GFED4s

472

estimates for 2024 are as of October 17.

473

474

S.3 Methodology Ocean CO2 sink

475

S.3.1 Observation

based estimates

476

We primarily use the observational constraints assessed by IPCC of a mean ocean CO

sink of 2.2 ± 0.7 GtC yr

477

for the 1990s (90% confidence interval; Ciais et al., 2013) to verify that the GOBMs provide a realistic

478

assessment of S

OCEAN

. This is based on indirect observations with seven different methodologies and their

479

uncertainties, and further using three of these methods that are deemed most reliable for the assessment of this

480

quantity (Denman et al., 2007; Ciais et al., 2013). The observation

based estimates use the ocean/land CO

sink

481

partitioning from observed atmospheric CO

and O

concentration trends (Manning and Keeling, 2006;

482

Keeling and Manning, 2014), an oceanic inversion method constrained by ocean biogeochemistry data

483

(Mikaloff Fletcher et al., 2006), and a method based on penetration time scale for chlorofluorocarbons (McNeil

484

et al., 2003). The IPCC estimate of 2.2 GtC yr

for the 1990s is consistent with a range of methods

485

(Wanninkhof et al., 2013). We refrain from using the IPCC estimates for the 2000s (2.3 ± 0.7 GtC yr

), and the

486

period 2002

2011 (2.4 ± 0.7 GtC yr

, Ciais et al., 2013) as these are based on trends derived mainly from

487

models and one data

product (Ciais et al., 2013). Additional constraints summarised in AR6 (Canadell et al.,

488

2021) are the interior ocean anthropogenic carbon change (Gruber et al., 2019) and ocean sink estimate from

489

atmospheric CO

and O

(Tohjima et al., 2019) which are used for model evaluation and discussion,

490

respectively.

491

We also use nine estimates of the ocean CO

sink and its variability based on surface ocean

maps obtained

492

by the interpolation of surface ocean

measurements. Seven of the methods cover a period from 1990

493

onwards due to severe restriction in data availability prior to 1990 (Figure 11), whereas two span the time period

494

from 1957 and 1959 onwards. These estimates differ in many respects: they use different maps of surface

495

different atmospheric CO

concentrations, wind products and different gas

exchange formulations as specified

496

in Table S3. We refer to them as

products. The measurements underlying the surface

maps are from

497

the Surface Ocean CO

Atlas version 2024 (SOCAT v2024; Bakker et al., 2024), which is an update of version

498

3 (Bakker et al., 2016) and the subsequent annual updates used in previous versions of the global carbon budget.

499

SOCAT v2024 has an additional 3.0 million

measurements with an estimated accuracy of better than 5

500

μatm relative to v2023. Of these, 2 million are from 2023 in a total of 210 data sets (Table S7), while the largest

501

addition from earlier years is from 2022 with 64 data sets new to SOCAT. For the 2023 data, there are a total of

502

178 data sets with measurements in the Northern hemisphere, while there are only 52 with data from the

503

Southern hemisphere. For the Southern Ocean, there are only 11 data sets from 2023 in the subpolar zone and

504

further south (defined as south of 45°S), and only one from Austral winter (June

August). The coverage of

505

SOCAT observations in 2023 is only about 50% of that in 2016 (Fig. 11), with large reductions in sampling in

506

both the Northern (from 391 to 178 data sets) as well as Southern hemisphere (from 109 to 52 data sets). This

507

reduction cannot be explained only in terms of lags in data submission. The quality control criteria used for

508

SOCATv2024 are described in Lauvset et al. (2018).

509

. Each of the data

based estimates uses a different method to map the SOCAT v2024 data to the global ocean.

510

The methods include a data

driven diagnostic method combined with a multi linear regression approach to

511

extend back to 1957 (Rödenbeck et al., 2022; referred to here as Jena

MLS), four neural network models

512

(Landschützer et al., 2014; referred to as VLIZ

SOMFFN; Chau et al., 2022; Copernicus Marine Environment

513

Monitoring Service, referred to here as CMEMS

LSCE

FFNN; Zeng et al., 2022; referred to as NIES

ML3;

514

Gregor et al. 2019, referred to as CSIR

ML6), one cluster regression approach (Gregor et al., 2024; referred to

515

OceanSODA

ETHZv2

), a multi

linear regression method (Iida et al., 2021; referred to as JMA

MLR), and one

516

method that relates the

misfit between GOBMs and SOCAT to environmental predictors using the extreme

517

gradient boosting method extending back to 1959 (Gloege et al., 2022).. The ensemble mean of the

based

518

flux estimates is calculated from these eight mapping methods. Further, we show the flux estimate of the

UExP

519

FNN

U method (Watson et al., 2020; Ford et al., accepted)

who also use a neural network model to map

520

data to the globe, but resulting in a substantially larger ocean sink estimate, owing to a number of adjustments

521

they applied to the surface ocean

data. Concretely, these authors adjusted the SOCAT

downward to

522

account for differences in temperature between the depth of the ship intake and the relevant depth right near the

523

surface, and included a further adjustment to account for the cool surface skin temperature effect. In

524

Friedlingstein et al. 2023, the UExP

FNN

U product correction was applied illustrating that this temperature

525

adjustment leads to an upward correction of the ocean carbon sink, up to 0.9 GtC yr

, that, if correct, should be

526

applied to all

based flux estimates. This year, the updated UExP

FFN

U method applies a smaller

527

adjustment as proposed by Dong et al. (2022), who illustrate a smaller correction effect of 0.6 GtC yr

. The

528

impact of the cool skin effect on air

sea CO

flux is based on established understanding of temperature gradients

529

(as discussed by Goddijn

Murphy et al., 2015 and Woolf et al., 2016), and laboratory observations (Jähne and

530

Haussecker, 1998; Jähne, 2019), but in situ field observational evidence is lacking (Dong et al., 2022). The

531

UExP

FNN

U method is thus, similar to the

UExP

FNN

U flux estimate in previous editions, not included in the

532

ensemble mean of the

based flux estimates. This choice will be re

evaluated in upcoming budgets based

533

on further lines of evidence.

534

Typically,

products do not cover the entire ocean due to missing coastal oceans and sea ice cover. The

535

flux from each

based product is already at or above 99% coverage (either due to complete coverage

536

or a posteriori filling) of the ice

free ocean surface area in several products this year (

UExP

FNN

JMA

MLR,

537

VLIZ

SOMFFN, Jena

MLS,

OceanSODA

ETHZv2

), . The products that remained below 99% coverage of the

538

ice

free ocean (CMEMS

LSCE

FFNN, NIES

ML3,

UExP

FNN

U, CSIR

ML6

) were scaled by the following

539

procedure:

540

Since v2022 of the GCB we now scale fluxes globally and regionally (North, Tropics, South) to match the ice

541

free area (using the HadISST sea surface temperature and sea ice cover; Rayner et al., 2003):

542

퐹퐶

푂

234

678&3#

(

'()

)

+),'-.

⋅

퐹퐶

푂

234:;'

543

In the equation,

represents area, (1

–

ice) represents the ice free ocean, A

FCO2

region

represents the coverage of

544

the

product for a region, and FCO

region

is the integrated flux for a region.

545

We further use results from two diagnostic ocean models, Khatiwala et al. (2013) and DeVries (2014), to

546

estimate the anthropogenic carbon accumulated in the ocean prior to 1959. The two approaches assume constant

547

ocean circulation and biological fluxes, with S

OCEAN

estimated as a response in the change in atmospheric CO

548

concentration calibrated to observations. The uncertainty in cumulative uptake of ±20 GtC (converted to ±1σ) is

549

taken directly from the IPCC’s review of the literature (Rhein et al., 2013), or about ±30% for the annual values

550

(Khatiwala et al., 2009).

551

552

S.3.2 Global Ocean Biogeochemistry Models (GOBMs)

553

The ocean CO

sink for 1959

2023 is estimated using ten GOBMs (Table S2). The GOBMs represent the

554

physical, chemical, and biological processes that influence the surface ocean concentration of CO

and thus the

555

air

sea CO

flux. The GOBMs are forced by meteorological reanalysis and atmospheric CO

concentration data

556

available for the entire time period. They mostly differ in the source of the atmospheric forcing data

557

(meteorological reanalysis), spin up strategies, and in their horizontal and vertical resolutions (Table S2). All

558

GOBMs except one (CESM

ETHZ) do not include the effects of anthropogenic changes in nutrient supply

559

(Duce et al., 2008). They also do not include the perturbation associated with changes in riverine organic carbon

560

(see Section 2.10 and Supplement S.6.3).

561

Four sets of simulations were performed with each of the GOBMs. Simulation A applied historical changes in

562

climate and atmospheric CO

concentration. Simulation B is a control simulation with constant atmospheric

563

forcing (normal year or repeated year forcing) and constant pre

industrial atmospheric CO

concentration.

564

Simulation C is forced with historical changes in atmospheric CO

concentration, but repeated year or normal

565

year atmospheric climate forcing. Simulation D is forced by historical changes in climate and constant pre

566

industrial atmospheric CO

concentration.

567

The atmospheric CO

forcing file was updated in GCB2024 to ensure consistency with the atmospheric CO

568

growth rate reported in the GCB. Since January 1980, we use the CO

global growth rate reported by

569

NOAA/GML (Lan et al., 2024). In the period March 1958

December 1979, we use bias

adjusted values of the

570

global growth rate based on measurements of atmospheric CO

made by the Scripps Institution of

571

Oceanography at the Mauna Loa Observatory, Hawaii (Keeling et al., 1976; full period of coverage 1758

2024).

572

Bias adjustment of the Scripps data was performed in three sequential stages as follows:

573

●

First, to correct for differences in the mean atmospheric concentration of CO

at Mauna Loa versus the

574

globally averaged value, a constant of

0.231 ppm was added to all Scripps data to improve alignment

575

of the “CO

[trend]” values from the Scripps data with the “CO

[trend]” values from the global NOAA

576

data. The value of

0.231 ppm is the mean offset of “CO

[trend]” at Mauna Loa from the global

577

“CO

[trend]” value during 1980

2000.

578

●

Second, to correct for differences in the seasonality of atmospheric CO

concentrations at Mauna Loa

579

versus globally, we shifted monthly anomalies between CO

concentration data and “trend” values

580

backward in time by one month in the Scripps data. This specifically corrects for the fact that

581

peaks/troughs in the climatology of "CO

[monthly_observation]

[trend]” at Mauna Loa occur 1

582

month earlier than peaks/troughs in the climatology of "CO

[monthly_observation]

[trend]” in the

583

global data from NOAA. A one

month shift to the Scripps data was found to optimally align the

584

climatologies of "CO

[monthly_observation]

[trend]” in the Scripps and global data.

585

●

Third, to correct for the greater amplitude of seasonal anomalies at Mauna Loa from Scripps than the

586

global data from NOAA, we apply a monthly multiplier that dampens the magnitude of monthly

587

anomalies from “trend” values in the Scripps data. The monthly multiplier reduces values of

588

"CO

[monthly_observation]

[trend]” in the Scripps data to more closely match values of

589

"CO

[monthly_observation]

[trend]” in the NOAA global data.

590

591

For the period Jan 1750 to February 1958, we use bias

adjusted values of the global growth rate based on

592

measurements of atmospheric CO

from air trapped in ice at Law Dome (Joos and Spahni, 2008; full period of

593

coverage 1750

2004). Bias adjustments were made to improve alignment with the post

1980 time series of data

594

from Scripps and NOAA, and were performed in two sequential stages as follows:

595

●

First, a constant of 0.973 was added to all data from Law Dome to improve alignment with the Scripps

596

data (which had already been bias

corrected as described above). The constant of 0.973 is the mean

597

offset of CO

annual values (annual mean in the case of the Scripps data) in the period 1958

1979.

598

●

Second, the climatology of "CO

[monthly_observation]

[trend]” from the period 1958

2000 was

599

superimposed on the data from Law Dome (note that the 1958

2000 data includes both Scripps and

600

NOAA data, combined as described above). To achieve this, a spline interpolation was fitted to

601

downscale annual observations from CO

concentration from Law Dome to monthly values of

602

“CO

[trend]” and the climatological seasonality of "CO

[monthly_observation]

[trend]” from

603

1958

2000) was then added to the interpolated values of “CO

[trend]”.

604

605

To derive S

OCEAN

from the model simulations, we subtracted the slope of a linear fit to the annual time series of

606

the control simulation B from the annual time series of simulation A. Assuming that drift and bias are the same

607

in simulations A and B, we thereby correct for any model drift. Further, this difference also removes the natural

608

steady state flux (assumed to be 0 GtC yr

globally without rivers), which is often a major source of biases.

609

Note, however, that Gürses et al. (2023) questioned the assumption of comparable bias and drift in simulations

610

A and B as they compared two versions of FESOM

REcoM, and found a very similar air

sea CO

flux in

611

simulation A despite a different bias as derived from simulation B. This approach works for all model set

ups,

612

including IPSL, where simulation B was forced with variable historical climate changes (looping over a 10

year

613

forcing). This approach assures that the interannual variability is not removed from IPSL simulation A.

614

The absolute correction for bias and drift per model in the 1990s varied between <0.01 GtC yr

and 0.31 GtC

615

, with five models having positive biases, four having negative biases and one model having essentially no

616

bias (NorESM). The MPI model uses riverine input and therefore simulates outgassing in simulation B. By

617

subtracting a linear fit of simulation B, also the ocean carbon sink of the MPI model follows the definition of

618

OCEAN

. This correction increases the model mean ocean carbon sink by 0.07 GtC yr

in the 1990s. The ocean

619

models cover 99% to 101% of the total ocean area, so that area

scaling is not necessary.

620

621

S.3.3 GOBM evaluation

622

The ocean CO

sink for all GOBMs and the ensemble mean falls within 90% confidence of the observed range,

623

or 1.5 to 2.9 GtC yr

for the 1990s (Ciais et al., 2013) before and after applying adjustments. The GOBMs and

624

products have been further evaluated using the fugacity of sea surface CO

(

) from the SOCAT v2024

625

database (Bakker et al., 2016, 2024). We focused this evaluation on the root mean squared error (RMSE)

626

between observed and modelled

and on a measure of the amplitude of the interannual variability of the

627

flux (modified after Rödenbeck et al., 2015). The RMSE is calculated from detrended, annually and regionally

628

averaged time series of

calculated from GOBMs and

products subsampled to SOCAT sampling

629

points to measure the misfit between large

scale signals (Hauck et al., 2020). To this end, we apply the

630

following steps: (i) subsample data points for where there are observations (GOBMs/

products as well as

631

SOCAT), (ii) average spatially, (iii) calculate annual mean, (iv) detrend both time

series (GOBMs/

632

products as well as SOCAT), (v) calculate RMSE. We use a mask based on the minimum area coverage of the

633

products. This ensures a fair comparison over equal areas. The amplitude of the S

OCEAN

interannual

634

variability (A

IAV) is calculated as the temporal standard deviation of the detrended annual CO

flux time series

635

after area

scaling (Rödenbeck et al., 2015, Hauck et al., 2020). These metrics are chosen because RMSE is the

636

most direct measure of data

model mismatch and the A

IAV is a direct measure of the variability of S

OCEAN

637

interannual timescales. We apply these metrics globally and by latitude bands. Results are shown in Figure S2

638

and discussed in Section 3.6.5.

639

640

In addition to the interior ocean anthropogenic carbon accumulation (Section 3.6.5) and SOCAT

, we

641

evaluate the models with process

based metrics that were previously related to ocean carbon uptake. These are

642

the Atlantic Meridional Overturning Circulation (Goris et al., 2018, Terhaar et al., 2022, Terhaar et al., in

643

review), the Southern Ocean sea surface salinity (Terhaar et al., 2021, 2022, 2024, Hauck et al., 2023b), the

644

Southern Ocean stratification index (Bourgeois et al., 2022) and the surface ocean Revelle factor (Terhaar et al.,

645

2022, 2024).

646

647

We follow the methodology of previous studies wherever possible, particularly the RECCAP model evaluation

648

chapter (Terhaar et al.,2024). The Atlantic Meridional Overturning Circulation from the GOBMs is here defined

649

as the maximum of the Atlantic meridional overturning streamfunction at 26°N. This is compared to data from

650