Country-level estimates of gross and net carbon fluxes from land use, land-use change and forestry

The reduction of CO₂ emissions and the enhancement of CO₂ removals related to land use are considered essential for future pathways towards net-zero emissions and mitigating climate change. With the growing pressure under global climate treaties, country-level land-use CO₂ flux data are becoming increasingly important. So far, country-level estimates are mainly available through official country reports, such as the greenhouse gas inventories reported to the United Nations Framework Convention on Climate Change (UNFCCC). Recently, different modelling approaches, namely dynamic global vegetation models (DGVMs) and bookkeeping models, have moved to higher spatial resolutions, which makes it possible to obtain model-based country-level estimates that are globally consistent in their methodology. To progress towards a largely independent assessment of country reports using models, we analyse the robustness of country-level CO₂ flux estimates from different modelling approaches in the period 1950–2021 and compare them with estimates from country reports.

Our results highlight the general ability of modelling approaches to estimate land-use CO₂ fluxes at the country level and at higher spatial resolution. Modelled land-use CO₂ flux estimates generally agree well, but the investigation of multiple DGVMs and bookkeeping models reveals that the robustness of their estimates strongly varies across countries, and substantial uncertainties remain, even for top emitters. Similarly, modelled land-use CO₂ flux estimates and country-report-based estimates agree reasonably well in many countries once their differing definitions are accounted for, although differences remain in some other countries. A separate analysis of CO₂ emissions and removals from land use using bookkeeping models also shows that historical peaks in net fluxes stem from emission peaks in most countries, whereas the long-term trends are more connected to removal dynamics. The ratio of the net flux to the sum of CO₂ emissions and removals from land use (the net-to-gross flux ratio) underlines the spatio-temporal heterogeneity in the drivers of net land-use CO₂ flux trends. In many tropical regions, net-to-gross flux ratios of about 50 % are due to much larger emissions than removals; in many temperate countries, ratios close to zero show that emissions and removals largely offset each other. Considering only the net flux thus potentially masks large emissions and removals and the different timescales upon which they act, particularly if averaged over countries or larger regions, highlighting the need for future studies to focus more on the gross fluxes.

Data from this study are openly available via the Zenodo portal at https://doi.org/10.5281/zenodo.8144174 (Obermeier et al., 2023).

1 Introduction

The carbon dioxide (CO₂) flux from land use, land-use change and forestry (LULUCF) is a key component of the global carbon cycle (Pongratz et al., 2021), and the net CO₂ flux from LULUCF (f_LULUCF) has contributed about 10 %–15 % to the total anthropogenic CO₂ emissions in recent decades (Friedlingstein et al., 2022a). Under the efforts towards achieving the long-term temperature targets formulated in the Paris Agreement, the importance of f_LULUCF is expected to increase in the future (Fuss et al., 2018), potentially contributing ∼ 30 % of the global climate change mitigation needed in 2050 to reach the 1.5 ^∘C target via both emission reduction and CO₂ removal (Roe et al., 2019). Despite its outstanding importance, estimates of f_LULUCF remain subject to high uncertainty. For example, within the annual Global Carbon Budget (GCB) assessments of the Global Carbon Project (GCP), the net f_LULUCF has the highest relative uncertainty among all terms (∼ 60 % for the most recent decade; see Friedlingstein et al., 2022a for GCB2022).

At the global scale, such as in the assessments of the Intergovernmental Panel on Climate Change (IPCC) and the GCB, f_LULUCF is usually assessed by modelling approaches (Friedlingstein et al., 2022a), namely semi-empirical bookkeeping models (BKs) and dynamic global vegetation models (DGVMs), and more recently by merging bottom-up inventory estimates built up from country-level information (Grassi et al., 2018, 2023a). Global modelling approaches have the advantage of being globally consistent, thus enabling, for example, the analysis of the global carbon cycle. Given recent improvements in global modelling approaches, such as better representation of processes related to land management, vegetation physiology, and soil biogeochemistry as well as increased spatial resolutions (Arneth et al., 2017; Blyth et al., 2021; Pongratz et al., 2021), several analyses of f_LULUCF estimates using multiple approaches, including models, at country (e.g. Federici et al., 2017; Rosan et al., 2021; Schwingshackl et al., 2022; Grassi et al., 2023a) and regional scales (e.g. Bastos et al., 2020a; Petrescu et al., 2021; Tubiello et al., 2021; Nabuurs et al., 2022) have recently been conducted. Yet, their results might be less robust and consistent at finer spatial scales. Bottom-up approaches, such as national greenhouse gas inventories (NGHGIs), quantify f_LULUCF based on inventory data that countries regularly submit to the United Nations Framework Convention on Climate Change (UNFCCC). Even though NGHGIs should fulfil best-practice guidance from the IPCC (IPCC, 2006), the quality, methodological complexity and provision of data used by the reporting countries vary (Grassi et al., 2021, 2022). Additionally, the Food and Agriculture Organization of the United Nations (FAO) disseminates independent bottom up estimates of LULUCF emissions and removals, based on forest area and carbon stock information provided every 5 years by the countries to the FAO Global Forest Resources Assessment (FRA), complemented by geospatial information on biomass fires and peatland degradation (Tubiello et al., 2022).

Land-based climate change mitigation can be achieved via two levers, namely by decreasing gross emissions and by increasing gross removals (the latter is often referred to as negative emissions; e.g. Fuss et al., 2018). These gross fluxes, which add up to the net f_LULUCF, act on different timescales and are in reality often linked by common land-use practices (note that throughout the paper the term “CO₂ fluxes” generally refers to anthropogenic fluxes from land use). For example, if fewer wood products are harvested in forestry, this quickly leads to lower emissions but also to lower forest regrowth and thus lower removals, balancing towards net-zero fluxes in the longer-term (Gasser et al., 2022). In this study, we define gross fluxes as the sum of all fluxes related to certain LULUCF practices that typically lead to emissions or removals, respectively, in line with the Global Carbon Project. Gross emissions from LULUCF (CO₂ fluxes from the biosphere to the atmosphere, i.e. C source) are mainly related to cropland or pasture expansion resulting in the destruction of natural ecosystems by deforestation, forest and peatland degradation, as well as land-use practices, such as biomass burning or forest management causing the decay of harvested wood products (HWPs; Friedlingstein et al., 2022a). Gross removals from LULUCF (CO₂ fluxes from the atmosphere to the biosphere, i.e. C sinks) are associated with land-use changes such as afforestation and reforestation including the regrowth of secondary forest after agricultural abandonment, as well as land-use influences such as forestry cycles and restoration of other (non-forest) ecosystems. It is noteworthy that the CO₂ removal potential of LULUCF is deemed the most suitable and easily scalable option for negative emissions, potentially providing 25 % of net greenhouse gas emissions reductions by 2030, primarily via afforestation, reforestation and management of existing forests (Griscom et al., 2017; Gidden et al., 2022; Smith et al., 2023). This is reflected by the fact that negative emissions from LULUCF are already widely included in the nationally determined contributions (NDCs) for climate change mitigation (Grassi et al., 2017; Smith et al., 2023). While this clearly highlights the importance of separately estimating gross emissions and gross removals, global assessments of f_LULUCF usually consider net f_LULUCF only, and, as stated by Houghton (2020), little attention has been paid to its components. While some studies explicitly separated gross fluxes early on (Tubiello et al., 2015; Federici et al., 2015), GCB studies considered gross fluxes for the first time in 2020 (Friedlingstein et al., 2020), with the most recent GCB2022 estimating global anthropogenic gross emissions at 3.8 ± 0.7 GtC yr⁻¹ and gross removals at 2.6 ± 0.4 GtC yr⁻¹ for 2012–2021, thus being 2–4 times larger than the global net f_LULUCF of 1.2 ± 0.7 GtC yr⁻¹ in this period (Friedlingstein et al., 2022a).

The aim of this study is to provide a comprehensive national and regional comparison that integrates the different approaches and definitions around the world and complements previous studies on selective countries or regions. We investigate the robustness of net f_LULUCF estimates from models, namely from BKs and DGVMs, we provide estimates at country and regional level, and, additionally, we compare net estimates from models and inventories. To identify the drivers of the spatio-temporally varying net f_LULUCF estimations from these assessments, we present and discuss country-level gross f_LULUCF estimates as well as the net-to-gross flux ratio from BKs in addition to net flux estimates. This allows us a more nuanced analysis and identification of the land component processes that are relevant in specific regions and provides a quantification of land-based climate change mitigation potentials distinctly for the two levers “reducing emissions” and “increasing removals”. The need for such an assessment is underlined by the increasing political relevance of LULUCF fluxes for countries' emissions reduction pledges, for example, in support of the ongoing Global Stocktake and the Glasgow Leaders' Declaration on Forests and Land Use of the 26th UN Climate Change Conference of the Parties (COP26) in Glasgow.

2 Overview of different f_LULUCF assessment methods

We use data from BKs and DGVMs, NGHGIs, and FAOSTAT (the Statistical Division of the Food and Agricultural Organization) for our analysis. The approaches are briefly described below and explained in more detail in Appendix A.

2.1 Bookkeeping models (BKs)

BKs are explicitly designed to estimate LULUCF fluxes following land management and land-use changes by tracking the carbon content in soil, vegetation and product pools, i.e. stocks and fluxes of carbon in the land biosphere and between land and atmosphere. They combine spatial information on land-use activities with observation-based carbon stock densities and specific response curves of soil and vegetation carbon for each land-use conversion type (for more details refer to Sect. A1 and Pongratz et al., 2014). Using separate response curves for carbon release (emissions) and carbon uptake (removals), BKs model both gross fluxes explicitly, and net f_LULUCF is derived as the sum of the two gross fluxes. BKs do not model the fluxes from peatland fire and drainage, but these emissions are added on top of the BK estimates in this study according to the approaches used in the GCBs; compare Sect. A1. We use simulations by three BKs that were conducted for the GCB2022 (Friedlingstein et al., 2022a), namely (1) the “bookkeeping of land use emissions” model (hereafter BLUE22; Hansis et al., 2015), (2) the “Houghton and Nassikas” model (hereafter H&N22; Houghton and Nassikas, 2017) and (3) the compact Earth system model “OSCAR” (hereafter OSCAR22; Gasser et al., 2020).

2.2 Dynamic global vegetation models (DGVMs)

DGVMs are process-based models used to simulate the interaction between land surface and vegetation processes with the atmosphere. They approximate net f_LULUCF as the difference in net biome productivity (NBP) of a simulation including LULUCF (similarly to the BKs, using spatial information on land-use activities) and a simulation excluding LULUCF (the latter using a time-invariant pre-industrial land-use map; for more details refer to Sect. A2 and Obermeier et al., 2021). DGVM simulations using transient (observed) environmental forcing are operationally available and can be used to estimate the impact of climate variability or long-term environmental changes on f_LULUCF. In addition, DGVM simulations can be forced with constant environmental data prescribing, for instance, pre-industrial or present-day environmental conditions. Simulations using the latter setup most closely resemble conditions that occurred during the time when observed carbon densities (as used by BKs or inventories) were measured and are therefore the best DGVM setup to compare with BKs or inventories. Here we use nine DGVMs that provided simulations with present-day as well as transient environmental forcing within the project “Trends and drivers of the regional-scale emissions and removals of carbon dioxide” (TRENDY; Le Quéré et al., 2014; Sitch et al., 2015).

2.3 National greenhouse gas inventories (NGHGIs)

NGHGIs are the official country reports to UNFCCC that include estimates of greenhouse gas emissions and removals from LULUCF. They widely rely on empirical emission factors in combination with country-level data on land-use activities to estimate f_LULUCF (for more details refer to Sect. A3 and Grassi et al., 2022). The methods and report details strongly vary between and among non-Annex I and Annex I countries. Non-Annex I countries frequently use default IPCC emission factors (IPCC, 2006) and often only report fluxes from deforestation, while estimates from Annex I countries are often based on country-specific statistical or process-based models for all IPCC land-use categories (forest land, cropland, grassland, wetlands, settlements and other land). For the analysis presented here, we use the country database (DB) compiled by Grassi et al. (2022), as updated in Grassi et al. (2023a) (hereafter referred to as NGHGI DB), which includes GHG data from all available country reports submitted to UNFCCC, gap-filled when necessary to allow a complete time series from 2000 until 2020. According to UNFCCC guidelines, country reports should encompass all LULUCF fluxes from any areas considered managed and for which IPCC provide methodological guidance. In most cases, NGHGIs include natural and indirect anthropogenic fluxes (e.g. from CO₂ fertilization; Grassi et al., 2018; Schwingshackl et al., 2022; IPCC, 2010). To improve comparability of the NGHGI DB data with modelled estimates, we adjust the NGHGI DB data to better match the processes and definitions of the modelled estimates (the basis in this study) by correcting for this so-called managed land issue (hereafter adjusted NGHGI DB). Following the approach described in Grassi et al. (2023a), we therefore subtract from the NGHGI DB estimates those fluxes resulting from natural and indirect effects (i.e. human-induced environmental change) from managed land. These effects are estimated by the ensemble mean of 16 transient DGVM simulations without land-use changes from TRENDYv11 (corresponding to the “natural terrestrial sink” in recent GCBs of the GCP), except for Brazil and Canada filtered with an intact/non-intact forest mask (Potapov et al., 2017). For Brazil and Canada, this approach uses the national gridded maps on managed and unmanaged forests used in the respective NGHGIs (Brazil, 2020; Canada, 2021).

2.4 Statistical Division of the Food and Agricultural Organization (FAOSTAT)

FAOSTAT f_LULUCF estimates resemble the bottom-up approach of the UNFCCC data in the way that they are based on consistent underlying activity data – at grid cell or at country level – in combination with emission factors. FAOSTAT f_LULUCF data are estimated by applying IPCC Guidelines (IPCC, 2003, 2006) to activity data generated either through official country reporting processes or through geospatial data analysed by FAO (for more details refer to Sect. A3 and FAO, 2020; Tubiello et al., 2021). FAOSTAT f_LULUCF cover carbon emissions and carbon removals in forests and from deforestation, derived from national carbon stock change statistics and IPCC emission factors, as well as emissions from peatland fires and peat drainage, the latter two obtained from satellite imagery in combination with IPCC emission factors (Conchedda and Tubiello, 2020; Tubiello et al., 2021; Prosperi et al., 2020). The IPCC (2006, 2019) recognizes explicitly that both FAO activity data and FAO emissions estimates provide a valuable set of reference data that can be used for validation, quality control and quality assurance of the data submitted through NGHGIs.

2.5 Main differences between the approaches

Although the different approaches summarized above (and described in detail in Appendix A) all aim at quantifying CO₂ fluxes from LULUCF, they differ substantially. Some of the key differences between the approaches are briefly described below; further details are provided in the results in Sect. 4.2 and 4.3 and in the Appendix.

Uncertainties in f_LULUCF estimations from modelling approaches mainly arise at high spatial resolutions (Kondo et al., 2022), from differences in underlying land-use and land-cover information (Gasser et al., 2020; Hartung et al., 2021; Ganzenmüller et al., 2022), missing observational constraints (Goll et al., 2015; Li et al., 2017), differences in process complexity and in the degree of implementation of LULUCF practices (Arneth et al., 2017; Hartung et al., 2021; Fisher and Koven, 2020), inconsistencies in common terminology and definitions (Pongratz et al., 2014; Gasser and Ciais, 2013), and different model assumptions and setups (Obermeier et al., 2021; Bastos et al., 2021a).

Most of the investigated modelling approaches use the spatially explicit LUH2-GCB2022 dataset as LULUCF forcing (Friedlingstein et al., 2022a). The BK model H&N22 uses FAO activity data at country level, and OSCAR22 uses both LUH2-GCB2022 and FAO activity data. To analyse the impact of different land-use forcing data, we further use BLUE data from the GCB2019 (hereafter BLUE19). As the BLUE model code was not changed between the GCB2019 and GCB2022, this allows the direct impact of the LULUCF forcing data to be isolated, which for these GCBs was based on HYDE3.2 and HYDE3.3, respectively. The main innovations in HYDE3.3 were the provision of yearly output from 1950 onwards, the update of the onset of agriculture based on new radiocarbon data and archaeological expertise indicating more spatial heterogeneity, the use of the latest satellite data with increased spatial resolution on an annual basis from 1992–2018 from the European Space Agency (ESA) and MapBiomas data for Brazil for the period 1985–2020, and the inclusion of more sub-national cropland and pasture data.

BKs do not explicitly represent biogeochemical processes and do not directly use environmental forcing data and thus usually exclude effects from climate change and meteorological or climate variability. DGVMs, in contrast, incorporate biogeochemical processes implemented in their modelling scheme, and – additionally to land-use change data – they use environmental forcing data as input. Consequently, DGVM estimates under transient environmental conditions include the long-term response to changing environmental conditions (e.g. atmospheric CO₂ increase, nitrogen deposition and fertilizer applications) and effects of climate variability. Due to the setup to calculate LULUCF emissions from DGVMs, transient DGVM f_LULUCF estimates include the loss of additional sink capacity (LASC), representing carbon fluxes in response to environmental changes on managed land (typically croplands with low carbon sink capacity and fast turnover rates) as compared to potential natural vegetation (typically forests with large carbon sinks and low turnover rates), which leads to higher flux estimates compared to BKs towards the end of the simulated periods (Gasser and Ciais, 2013; Pongratz et al., 2014; Obermeier et al., 2021).

In contrast to modelling approaches that provide globally consistent historical f_LULUCF analysis over the entire simulated period (e.g. from 1700 onward), country reports only cover the most recent decade for which observations and statistics exist (starting around 1990–2000, depending on the dataset and country) and are restricted to the territories and components that countries report. In addition, f_LULUCF estimates from country reports rely on different definitions and assumptions than global models and therefore do not allow, for example, a realistic derivation of countries' climate change mitigation contributions that are consistent with the pathways to achieve the climate targets of the Paris Agreement as derived by integrated assessment models. One major difference is the definition of managed land, with NGHGIs having comparatively larger areas of managed land than models (Grassi et al., 2023a). Further, NGHGI estimates that rely on direct observations (e.g. national forest inventories) include both direct and most indirect anthropogenic effects on managed lands, which can only be separated based on modelling approaches (Grassi et al., 2018; Schwingshackl et al., 2022). Consequently, most NGHGI f_LULUCF estimates include large parts of the natural response to recent environmental change in managed lands (e.g. larger vegetation carbon sink due to so-called CO₂ fertilization) and, thus, estimate larger anthropogenic CO₂ removals (and widely lower net f_LULUCF) than global models (IPCC, 2006; Grassi et al., 2018; Schwingshackl et al., 2022). In order to make the NGHGI and BK land-use flux estimates comparable, here we translate NGHGI estimates by removing the fluxes that models attribute to the natural land sink (compare the adjusted NGHGI DB data described in Sect. A3).

Reported f_LULUCF estimates remain highly uncertain, particularly in many developing countries, as country reports to UNFCCC do not explicitly separate managed from unmanaged forest land or only report forest-related fluxes (Grassi et al., 2022). Other countries report fluxes from additional LULUCF practices, such as from natural peatland that is converted to agriculture or from land that is converted to human settlements (not included in modelled estimates). While human-induced degradation from logging and fires is often included in national reports, forest degradation fluxes are hardly existent in BKs.

The main differences between FAOSTAT data and the NGHGI database are explained by a more complete coverage in the NGHGI database than FAOSTAT, in particular the inclusion of non-biomass carbon pools and non-forest land uses (Grassi et al., 2022). Additional differences stem from (1) differences in activity data, for instance, different time series of the forest area which may be communicated independently to UNFCCC and FAO by different national agencies; (2) differences in emission factors, including carbon stock data and their sub-national resolution as well as their changes over time; and (3) differences in scale, especially considering that forest fluxes may be computed at grid cell or sub-national scale in many countries using remote sensing or information from national forest inventories. Additionally, the FAO area data do not distinguish between managed and unmanaged areas and, thus, do not in principle separate anthropogenic and non-anthropogenic drivers (Tubiello et al., 2021; Grassi et al., 2022). The nature of data reported to FAO indicates that in most cases the area considered in FAOSTAT estimates is similar to the one in the NGHGIs, but often the effects of environmental changes (i.e. natural and indirect anthropogenic effects) are not included as in the NGHGIs. For this reason, we did not correct FAOSTAT estimate data here.

3 Data processing for country-level and regional analysis

This study compiles f_LULUCF estimates from various modelling and country-report-based approaches, aggregated to the country and regional level. Country-level aggregation is provided for all 186 (out of 195) UNFCCC country parties that reported LULUCF fluxes, comprising > 99.6 % of global net f_LULUCF (as derived by the three BKs). Regional aggregation is based on the 18 land regions defined by the REgional Carbon Cycle Assessment and Processes Phase 2 (RECCAP2; part of European Space Agency Climate Change Initiative (ESA CCI); Ciais et al., 2022), as defined in Tian et al. (2019) and shown in Fig. A1.

Due to different spatial resolutions of the investigated datasets, preliminary processing steps were needed to obtain f_LULUCF at country level. Data from NGHGI DB, FAOSTAT, H&N22 and OSCAR22 are disseminated at country level already. For the NGHGI DB and FAOSTAT estimates, we converted CO₂ fluxes into carbon fluxes based on their molar mass fraction, i.e. dividing by $\mathrm{44}/\mathrm{12}$ (according to the IPCC (2006) approach). The gridded outputs of the DGVMs and of the BLUE model were aggregated to country-level f_LULUCF estimates based on a (modified) 0.25^∘ country mask from Columbia University – Center for International Earth Science Information Network (CIESIN) (2018). We remapped the mask of each country to the native grid of each DGVM by conservative remapping using Climate Data Operators (CDO), which yields the area fraction of every country in each DGVM grid cell. We calculate each country's f_LULUCF share by multiplying the total f_LULUCF in a grid cell by the share of each country of the total land fraction in that grid cell. The country-wide f_LULUCF is then obtained as a sum over all grid cells.

Similar to the country-level aggregation, the RECCAP2 region mask was remapped to the native grid of each DGVM and the BLUE model to obtain regional aggregation. Country-level data from FAOSTAT, H&N22, NGHGI DB, and OSCAR22 were aggregated to RECCAP2 regions by overlaying the CIESIN country map. Where country borders did not match RECCAP2 regions, namely where a country is split between two regions, those country's fluxes in FAOSTAT, H&N22, NGHGI DB, and OSCAR22 were attributed to the RECCAP2 containing the largest area fraction of the country.

4 Country-level f_LULUCF estimates from different approaches

In the main paper we exemplarily focus on the top eight countries with highest cumulative net LULUCF emissions in the period 1950–2021 based on BKs, as well as the United States with the highest cumulative removals in this period. In the following, we analyse and compare f_LULUCF estimates from BKs, DGVMs, and country-report-based approaches to assess the robustness of f_LULUCF estimates within and across the different datasets. Specifically, we identify countries in which the different estimates agree well and those where uncertainties in f_LULUCF estimates are high. We further quantify and assess the gross fluxes (i.e. gross emissions and gross removals) in each country and analyse their importance in comparison to net fluxes.

Appendix A contains a comprehensive compilation of the f_LULUCF estimates from BKs, DGVMs, and country-report-based approaches for the RECCAP2 regions (Figs. A4–A8) and for the investigated 186 countries (Figs. A9–A18). Summary statistics of the BK estimates for all 186 country aggregates and the EU27 + UK can be found in Tables A1–A3, and a dataset covering all aggregated country data for the period 1950–2021 for all datasets used is accessible under https://doi.org/10.5281/zenodo.8144174 (Obermeier et al., 2023).

4.1 Overview of selected country-level net f_LULUCF estimates from bookkeeping models

Country-level net f_LULUCF estimates from bookkeeping models strongly vary across countries. Net LULUCF emissions are highest (in descending order) in Brazil, Canada, China, DR Congo, India, Indonesia, Nigeria and Russia based on cumulative estimates between 1950 and 2021 (Fig. 1a). The eight countries with the highest cumulative net emissions from LULUCF are either located in carbon-rich, forested tropical regions (Brazil, DR Congo, Indonesia and Nigeria) and/or encompass vast territories (Brazil, Canada, China, India and Russia). These top eight emitters emitted more than ∼ 53 % of the total net LULUCF emissions in the period 1950–2021 and are thus of outstanding importance for climate change mitigation via LULUCF emission reductions. In particular in Indonesia, the estimated emissions per area are higher than in all other 185 countries studied, which illustrates an enormous pressure on the terrestrial carbon stocks in the tropics (Fig. A3). However, it should be noted that the resulting call for action is not limited to these top emitters, as large parts of national LULUCF emissions, particularly in the tropics, are caused by consumption elsewhere (Hong et al., 2022).

Figure 1 Net carbon fluxes from land use, land-use change and forestry (f_LULUCF) from three bookkeeping models (BKs; data from GCB2022 simulations). (a) Cumulative carbon fluxes over 1950–2021 and (b) average carbon fluxes in 2011–2021. The bars show the mean of the three BKs (filled bars) and minimum and maximum estimates (hatched bars). Numbers in parentheses show the multi-model average and standard deviation (in GtC in (a) and MtC yr⁻¹ in (b)). Colours indicate the absolute quantities, showing countries with net emissions in red and countries with net removals in green. All 186 country aggregates from this study are shown in decreasing order of their (a) cumulative and (b) most recent annual f_LULUCF. In each panel, the top 10 emitters and the five countries with the largest removals are labelled (and the countries from the main paper if not yet included). The dashed red lines show the percentiles of net carbon emissions for each panel when adding the countries in decreasing order.

During the second half of the 20th century, high net LULUCF hotspots became increasingly concentrated in the countries of the Global South. As stated in the GCB2022 (Friedlingstein et al., 2022a), more than 50 % of most recent net emissions from LULUCF occur in only three countries, namely Brazil, DR Congo and Indonesia – all located in the tropics (Fig. 1b). The pronounced land-use change impact in the tropics is also evident from the f_LULUCF estimates calculated per area and per capita, with the 10 largest net emitters all located in the tropics (Figs. A3 and A2). However, most of these emissions are embodied in trade and are caused by consumption in industrialized regions such as Europe, the United States and China (Hong et al., 2022). A trend towards fewer countries comprising larger shares of global net LULUCF emissions is found when comparing cumulative and most recent annual f_LULUCF estimates. In the period 1950–2021, the top 22 net emitters comprised more than 75 % (top 43 emitters more than 90 %) of cumulative net emissions from LULUCF, while in 2011–2021, only 15 countries make up 75 % (top 33 emitters more than 90 %) of the global net emissions. In China, the f_LULUCF estimates from the BKs show an remarkable turnaround, turning China from the third-highest cumulative emitter in 1950–2021 to the country with the third-highest net removals in 2011–2021.

Table 1 Statistics of the annually averaged (2011–2021) and cumulative (1950–2021) net carbon fluxes from land use, land-use change and forestry (f_LULUCF) for the nine countries analysed. The table indicates the mean f_LULUCF estimates, their standard deviation across the different model estimates (SD) and their relative uncertainty (SD divided by the absolute mean value).

The large number of net emitting countries and the fact that BKs estimate a global net source of carbon from LULUCF to the atmosphere are in stark contrast to the pledges to achieve the goals associated with the Paris Agreement and the reported global carbon removal from LULUCF when summed across all NGHGIs (Grassi et al., 2022). Of note, despite widely included net negative emissions from LULUCF in the NDCs, including the need for carbon dioxide removal (CDR) techniques, BKs estimate net removals from LULUCF only for very few countries (both cumulatively between 1950–2021 as well as more recently in 2011–2021). Substantial country-level net removals from LULUCF are only found in the United States; some European countries; and, in the most recent decade, in China. However, large relative uncertainties in BK estimates remain, within the cumulative net f_LULUCF estimates (globally ∼ 38 %, and much higher in specific countries, for example, in China ∼ 150 %, the United States ∼ 680 %; see Table 1) and recent annual f_LULUCF estimates (globally ∼ 31 %, in China ∼ 700 % and Russia ∼ 1000 %). In order to explain the uncertainties related to f_LULUCF estimates from BKs, the latter need to be compared to estimates from other approaches, and the underlying drivers for the differences in f_LULUCF estimates need to be investigated.

4.2 Comparing country-level net f_LULUCF estimates from BKs and DGVMs

To test the robustness and explain the large spread in modelled country-level f_LULUCF estimates, we compare and discuss the differences of the net estimates from the BK ensemble and the DGVM ensemble with respect to the characteristics of the individual modelling approaches, particularly regarding the underlying land-use forcing data and, for the DGVM ensemble, the environmental forcing data.

Figure 2 Time series of net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2021 derived by bookkeeping models (BKs) and TRENDYv8 simulations with dynamic global vegetation models (DGVMs; used in the GCB2019) under present-day climate forcing. The eight countries with highest cumulative emissions since 1950 are shown in decreasing order and the United States. The figure shows the mean and absolute range of three BKs (using the GCB2022 simulations, BKs 2022) and the median and interquartile range of the eight DGVMs. Additionally, estimates from BLUE simulations from the GCB2019 (BLUE19; blue solid) and from the GCB2022 (BLUE22; blue dashed) are shown to illustrate the impact of updates in the LUHv2 forcing data (difference shows in steel blue) and the differences between DGVMs and BLUE (TRENDYv8 present-day minus BLUE19; purple) to illustrate the relevance of different modelling approaches. Numbers in the top-right corner of each panel indicate the median of the cumulative f_LULUCF sums over 1950–2018 (in GtC) for TRENDYv8 (red), BLUE19 (blue), BLUE22 (blue), and the difference of TRENDYv8 present-day minus BLUE19 (purple) and of the LUHv2 forcing (BLUE22 minus BLUE19; steel blue); numbers in the bottom-left corner indicate coefficients of correlation squared for BLUE19 and TRENDYv8 (purple) and for BLUE19 and BLUE22 (steel blue). BLUE19 data are only available until 2019, and TRENDYv8 data are only available until 2018. Greenish background depicts negative f_LULUCF, which is carbon removal from the atmosphere.

Country-level net f_LULUCF estimates of BKs and DGVMs from 1950 onward agree generally well in the investigated countries (Fig. 2; refer to Fig. A4 for RECCAP2 regions and Figs. A9 and A10 for all 186 investigated countries).

In most countries, modelled estimates from BLUE19 and (present-day) DGVM simulations show consistent temporal evolutions and peaks in emissions, when using the same land-use forcing data (compare the “difference of the modelling approach”, derived as the present-day TRENDYv8 mean minus BLUE19 which share the same (LUHv2) land-use forcing data; purple line in Fig. 2). These consistent trends and peaks lead to comparably high correlation coefficients between the BLUE19 and present-day DGVM estimates, for example, in Brazil, China and Nigeria. However, the correlation coefficients for the estimates from the different modelling approaches are generally rather low, which is due to the interannual variability, which is captured by the DGVMs but usually not by the BKs. In line with this, the particularly low correlation coefficients between DGVMs and BKs in the United States and Russia result from high interannual variability in combination with a low signal from land-use changes. Yet, the f_LULUCF estimates also substantially differ in some of the countries with pronounced land-use changes, despite identical land-use forcing data. The most striking difference occurs in 1997 in Indonesia. The year 1997 was an El Niño year, causing high carbon emissions from organic soils in Indonesia, which are included in the BK estimates but not in the DGVM estimates used (in line with the GCP assessments).

BLUE19 generally yields higher net emission estimates compared to the present-day TRENDYv8 ensemble (except for Nigeria), which indicates a tendency towards higher estimates when using the BK approach. To investigate this further, we compare the estimates of all three BKs from 2022 (BKs 2022; note that BLUE model code was not changed between 2019 and 2022 versions). BLUE22 tends to estimate the highest emissions among the BKs in most countries and during most of the time, indicating that the BK mean might agree better with the mean of the present-day DGVM simulations compared to BLUE19. This can mainly be explained by higher differences in the carbon densities between natural and managed areas assumed in the BLUE model in conjunction with the fact the BLUE captures the full extent of (LUH2-based) gross transitions (compare Sect. 4.3 and description of the BKs in Sect. A1 and Bastos et al., 2021a).

The strong influence of the land-use forcing data is highlighted by the often differing trends and peaks in f_LULUCF from BLUE19 based on HYDE3.2 and BLUE22 based on HYDE3.3 (compare the “LUHv2 forcing difference”, derived as BLUE22 minus BLUE19, in Fig. 2) and the huge ranges in the BK estimates in many regions. Large LUHv2 forcing differences, as indicated by particularly low correlation coefficients, are found in Brazil, DR Congo and Nigeria, with substantially larger estimates using the more recent HYDE3.3 data. The particularly big differences for Brazil mainly result from an improved representation of deforestation patterns through the inclusion of MapBiomas land cover data in the newer HYDE version, and, for DR Congo, they result from the inclusion of revised data from the FAO (Friedlingstein et al., 2022b). In contrast, BLUE22 has lower f_LULUCF estimates than BLUE19 in Canada (mainly in the 1950s) and, in the most recent decade, in China, India, Indonesia and Nigeria. In particular for tropical countries, lower emissions may be related to decreased cropland expansion in the updated HYDE data (Friedlingstein et al., 2022b).

Emission peaks around 1960 occur in several regions (e.g. Canada, DR Congo and Russia), mainly in the 2022 estimates (BLUE22) but not in the 2019 estimates (BLUE19). Those peaks can be explained by an artefact in HYDE3.3, resulting from the merge of historical HYDE data (used up to 1960) with FAOSTAT data (used from 1961 onward) (Chini et al., 2021; Friedlingstein et al., 2022b). The 1960 peaks thus do not represent any real-world land-use change fluxes and should be corrected in future updates of HYDE (as has already been achieved for Brazil; compare Sect. A1 and Friedlingstein et al., 2022b; Rosan et al., 2021). For Brazil, an additional peak in 2004 corresponds to increased deforestation for cropland and pasture, followed by a slowdown in deforestation rates due to governmental regulations, which is only captured in HYDE3.3 through the inclusion of ESA CCI Land Cover data (Rosan et al., 2021).

Figure 3 Net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2021 derived from TRENDYv8 simulations with dynamic global vegetation models (DGVMs) under historical (transient) and fixed present-day environmental conditions (compare Appendix). The eight countries with the highest cumulative emissions since 1950 are shown in decreasing order and the United States. The median and interquartile range of eight DGVMs are shown. Lines (shading) indicate the median (interquartile range) of the eight DGVMs, for which data for both simulations are available. Lines in the lower panels indicate the impact of the environmental forcing (present-day minus transient; purple). Numbers in the top-right corner depict the multi-model median of the cumulative sums in the period 1950–2019 (left column; GtC) and the differences between the simulations (right column); numbers in the bottom-left corner indicate the coefficient of correlation squared for TRENDYv8 present-day and transient simulations (purple). The light-green background indicates negative f_LULUCF, that is net carbon removal from the atmosphere by LULUCF.

Additionally, the environmental forcing plays an important role in many regions, as can be assessed by comparing transient and present-day TRENDYv8 simulations (Fig. 3 – the “Present-day versus transient difference”; refer to Fig. A5 for RECCAP2 regions and Figs. A11 and A12 for all 186 investigated countries). Estimates of f_LULUCF under present-day environmental forcing tend to be higher compared to transient estimates in the earliest simulated periods, particularly in tropical regions with high LULUCF activity (e.g. in Brazil, China, DR Congo and India), which is also reflected in the cumulative emissions estimates (see numbers displayed in Fig. 3). This can be explained by higher carbon stocks under present-day environmental forcing compared to transient forcing (the multi-DGVM mean global vegetation carbon stock increased by ∼ 23 % from 664 to 815 PgC from 1800 until today), particularly in early simulation periods (when atmospheric CO₂ concentration was substantially lower than today; Obermeier et al., 2021). Differences between present-day and transient f_LULUCF estimates become smaller as the simulations progress, as transient environmental conditions approach present-day conditions and additionally accumulate the loss of additional sink capacity (Pongratz et al., 2014). It is noteworthy that towards the end of the simulated period, transient f_LULUCF estimates even tend to exceed present-day estimates, due to the steadily accumulating loss of additional sink capacity, which globally comprises ∼ 0.8 ± 0.3 GtC yr⁻¹ (∼ 40 %) of transient f_LULUCF in the period from 2009–2018 (Obermeier et al., 2021). In the temperate zone (e.g. in Russia and the United States), variations in f_LULUCF estimates due to environmental conditions rather depend on the inter-annual meteorological and climate variability. Here, carbon stocks have not been enhanced by higher CO₂ concentrations as homogeneously as in the tropics, and climate change has even led to decreased carbon stocks in some regions (Obermeier et al., 2021). What is striking here is the particularly low correlation coefficient in the United States. Similar to the difference between the modelling approaches, this results primarily from the combination of a low land-use signal with very high interannual variability in environmental conditions.

As described above (and in the Appendix), the multiple modelling approaches differ in their underlying assumptions, their implemented process complexity and parametrizations, and the input forcing data used. These differences partly lead to highly differing estimates, in particular at finer spatial scales, which decreases the accuracy of some country-level estimates from models. The differences in country-level net f_LULUCF estimates that are due to the modelling approach versus changes in land-use forcing (approximated by the LUHv2 forcing difference for the BK model BLUE) or environmental forcing (difference between present-day and transient DGVM simulations) are of a similar order of magnitude, yet with very different spatial patterns. Whether the modelling approach, land use or environmental forcing is the dominating factor depends on the specific country and partly also on the specific time period. The modelling approach had the highest influence on the cumulative estimates, for example, in Brazil, Canada, India, Indonesia and Russia, whereas the LUHv2 forcing difference in the BK model BLUE was higher in China, DR Congo and Nigeria. For the DGVMs, the land-use forcing data impacted cumulative f_LULUCF estimates widely more strongly compared to the environmental forcing (except in the United States), although of similar magnitude. It should be emphasized that the strong influence of environmental factors, which is reflected in the up- and downswings of the DGVM estimates, is likely to become more important under future climate conditions with more frequent and intense extreme environmental conditions and potentially decreasing LULUCF activities as set out by the Glasgow Leaders' Declaration on Forests and Land Use.

4.3 Net f_LULUCF from individual bookkeeping models at country-level and country-report-based estimates

Despite broad agreement in most country-level net f_LULUCF estimates when averaged across multiple models, individual BK model output differs strongly in some countries (compare BK range in Figs. 2 and 4), and, not surprisingly, country-report-based estimates mostly diverge even more (Fig. 4; refer to Fig. A6 for RECCAP2 regions and Figs. A13 and A14 for all 186 investigated countries).

Figure 4 Net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2019 derived from individual bookkeeping models (BKs) used in GCB2022 (BKs 2022) and country-report-based estimates (FAOSTAT and NGHGI DB). The use of dashed and solid lines indicates that BK estimates and country-report-based estimates are not directly comparable (see Sect. 2.5 and the Appendix). For better comparability, the adjusted NGHGI DB estimates, matching the f_LULUCF definition of the BKs, are additionally shown. The grey line (shading) depicts the median (range) of the three BKs. The light-green background indicates negative f_LULUCF that is net carbon removal from the atmosphere by LULUCF.

Differences in annual net f_LULUCF estimates across the three BKs (based on simulations with the most recent land-use forcing) are particularly high in Canada and the United States throughout the 20th century and before and after emission peaks (e.g. 1960 and 2011 in DR Congo, 1980s in China, and in the 1950s in India). As stated above, the high differences related to emission peaks are predominantly due to the use of different land-use forcing data, which is reflected by the fact that H&N22 (based on FAO data) does not capture any of these peaks, while BLUE22 and OSCAR22 models (both of which use LUH2 data) show very similar trends and peaks. In China, for example, the H&N22 land-use forcing data assume a steady increase in forest areas from 1950, while the LUH2 data show decreasing forest areas until 1990 and relatively stable forest areas thereafter (Yu et al., 2022). The estimates of OSCAR22 often lie close to the multi-model mean of the BKs, which can be explained by the fact that OSCAR22 uses both FAO and LUH2 forcing data to derive a best-guess estimate for f_LULUCF (Gasser et al., 2020). The huge uncertainties in Canada, China and the United States cannot fully be explained by the net f_LULUCF but are discussed in Sect. 4.4 in relation to the gross fluxes.

Beyond the effects of the land-use forcing data, BK model differences in net f_LULUCF can be explained by several individual model specifics. The upper limit of the BK range is predominately defined by BLUE22, particularly during phases of high BK uncertainty, while the lower limit is often defined by H&N22. This is mainly due to different assumed carbon densities, which, depending on the ecosystem, are particularly high in the BLUE model and low in the H&N model (compare Appendix and Bastos et al., 2021a, b). Moreover, the inclusion of sub-grid-scale transitions in BLUE22 (and OSCAR22) increases emission estimates compared to H&N22. Additionally, H&N22 assumes conversion of natural grasslands to pasture, while BLUE22 and OSCAR22 allocate pasture proportionally on all natural vegetation that exists in a grid cell (Hansis et al., 2015; Friedlingstein et al., 2022b), yielding lower f_LULUCF for H&N22. In addition, different turnover periods for the HWPs in the different BKs lead to varying BK estimates after significant changes in LULUCF practices. The exception of slightly higher emission estimates from H&N22 in Indonesia (mainly in the 1990s) and DR Congo (from 1970s until 2007) is due to much higher carbon removals due to afforestation and reforestation in BLUE22 and OSCAR22 (not shown) and faster increasing emissions from deforestation in the H&N22 model.

Country-report-based net f_LULUCF estimates are substantially lower than BK estimates in most of the investigated countries; in particular NGHGI DB has the lowest emission estimates in almost all investigated countries. Much of this discrepancy (globally adding up to ∼ 1.6 GtC yr⁻¹ for the period 2001–2020) can be explained by different definitions, particularly regarding two points (compare Appendix and Grassi et al., 2023a; Schwingshackl et al., 2022). (1) Natural and indirect human-induced fluxes are included, such as those resulting from increased forest regrowth due to higher atmospheric CO₂ concentration and N deposition, in many country reports. (2) The area assumed to be managed land is larger and is therefore included in the f_LULUCF assessment in the country reports compared to the BKs. This is confirmed by the fact that the adjusted NGHGI DB data, where the natural land sink on managed land is subtracted, agree better with the BK estimates than the NGHGI DB data for most countries.

However, for individual countries, other reasons, such as omitted fluxes due to incomplete reporting of land uses and carbon pools, also (partly) explain the differences between modelled and reported f_LULUCF estimates (Schwingshackl et al., 2022). In the United States, lower country-report-based estimates may additionally result from the inclusion of CO₂ removals from urban vegetation (Churkina et al., 2010), which is not considered in the BK estimates. In addition, despite the increased methodological complexity of the approaches being used by many developed countries, it was shown that reported f_LULUCF estimates for the United States still remain uncertain, mainly due to uncertainty in inputs, model parameters and plot-based sampling (e.g. McGlynn et al., 2022). The Canadian NGHGI report discounts fluxes (emissions and removals) on areas affected by wildfires and severe insect disturbances and reports them in a separate category (Kurz et al., 2018), which can explain the strongly increased difference for the adjusted NGHGI DB estimate (Grassi et al., 2018; Schwingshackl et al., 2022). In China, the large gap between country-report-based and modelled estimates might be explained by high carbon removals from afforestation and ecological restoration projects (Yang et al., 2022; Jin et al., 2020; Yu et al., 2022) considered in the country reports but not fully included in the land-use data used by the models (compare Sect. 4.4 and Yu et al., 2022).

The generally good agreement in the f_LULUCF estimates from BKs and inventories in Indonesia is due to the dominance of the added peat data, with a pronounced interannual variability controlled mostly by El Niño–Southern Oscillation patterns on top of land use (Federici et al., 2017).

LULUCF flux estimates from FAOSTAT are largely in agreement with BK estimates, with some important differences, notably lower estimates in China and Russia and higher estimates in DR Congo and Canada. FAOSTAT estimates are generally higher than the NGHGI DB estimates except for the period from 2000–2019 in Russia and Nigeria and the emission peak in Brazil in 2004. As for the NGHGI DB estimates, lower FAOSTAT estimates compared to BK estimates are likely due to the inclusion of all fluxes on managed land (compare adjusted NGHGI DB data). Differences between FAOSTAT and NGHGI DB estimates can be explained by the generally more complete coverage of carbon fluxes in the latter and differing approaches to estimate forest fluxes, where FAO applies a carbon stock change approach based on observed forest data from FRA, while NGHGI reports are based on the use of a simple carbon stock change approach or a gain loss approach by the scaling up of forest growth rates based on IPCC default factors to forest land estimates (refer to Appendix and Grassi et al., 2022; Tubiello et al., 2021, for a detailed description of the differences between UNFCCC country, NGHGI DB and FAOSTAT estimates).

Additionally, the underlying data on forest land differ, with the NGHGI DB database reporting much greater forest areas and forest carbon removals (in particular for non-Annex I countries). Moreover, NGHGIs of Annex I and the largest non-Annex I countries also include non-biomass carbon pools and non-forest land uses, while, except for organic soils, FAOSTAT only includes above- and belowground biomass pools (Federici et al., 2015; Tubiello et al., 2021; Grassi et al., 2022). The higher estimates of NGHGI DB compared to FAOSTAT in Brazil are caused by larger deforestation and afforestation areas in the Brazilian report to UNFCCC compared to FRA and the fact that it considers gross deforestation and afforestation, while FAOSTAT reports net deforestation and afforestation directly (Federici et al., 2017; Rosan et al., 2021; Schwingshackl et al., 2022) (which might increase emissions in the NGHGI due to the asymmetry in instantaneously occurring gross emissions versus slowly increasing gross removals over the long term).

4.4 Gross f_LULUCF from bookkeeping models at country level

To get more insights into the underlying drivers of country-level net LULUCF estimates, we split them into gross fluxes, namely gross emissions (or “sources”) and gross removals (or “sinks”). As stated in the Introduction, here we define these gross fluxes as the sum of all fluxes related to those LULUCF practices that typically lead to emissions or removals, respectively. Gross emissions are mainly caused by deforestation, peatland degradation, biomass burning, the decay of HWPs and biomass left on site after harvest (Friedlingstein et al., 2022a). Gross removals are mainly associated with afforestation and reforestation including forest regrowth after agricultural abandonment, as well as forestry cycles and the restoration of other (non-forest) ecosystems. In the following, we present and discuss gross f_LULUCF derived from the three BKs for the nine selected countries (Fig. 5; refer to Fig. A7 for RECCAP2 regions and Figs. A15 and A16 for all 186 investigated countries).

Figure 5 Gross and net fluxes from land use, land-use change and forestry (f_LULUCF) from 1950 until 2021, as derived by three bookkeeping models (BKs; used in the GCB2022). Eight countries with the highest cumulative emissions since 1950 are shown in decreasing order and the United States. Lines (shaded area) depict the mean (range) of the BK estimates.

Similarly to the net fluxes, gross fluxes modelled by the three BKs show widely similar trends and agree on the timing of emission peaks. Peaks in net emissions are predominantly due to peaks in gross emissions, while the time series of gross removals is much smoother, consistent with the slower pace of vegetation regrowth. Consequently, the short-term evolution of net f_LULUCF is much more influenced by the dynamics of highly fluctuating gross emissions than by the dynamics of rather slowly changing gross removals. However, the decreasing trends in net f_LULUCF estimates across many regions (globally from 1.6 ± 0.7 GtC yr⁻¹ in the 1960s to 1.1 ± 0.7 GtC yr⁻¹ in the period from 2011–2020) mainly relate to a steady increase in the gross removals (globally from −1.9 ± 0.4 to −2.7 ± 0.4 GtC yr⁻¹ in the same period), which exceeded the increase in gross emissions (globally from 3.4 ± 0.9 to 3.8 ± 0.6 GtC yr⁻¹ in the same period; Friedlingstein et al., 2022b). This smooth evolution towards increased removals results from increased carbon sequestration on previously managed land, mainly due to forest regrowth and soil recovery. Environmental changes, such as those resulting from CO₂ fertilization, are not modelled by BKs (except the changing carbon densities in OSCAR).

In most regions, uncertainties in net f_LULUCF are due to uncertainties in gross emissions rather than uncertainties in gross removals. The largest differences, and thus pronounced uncertainties, in gross emissions from BKs are found for India and Canada as well as in the years before and after most emission peaks in many regions. Uncertainties related to the peaks in the 1960s mainly stem from the merging of two different datasets. In China, the pronounced peak in the 1980s is caused by spurious signals in the LUHv2 data, inherited from an abrupt cropland increment in the FAO data (Yu et al., 2022). Because cropland area is quantified relative to forest proportions, an increasing cropland area causes decreasing forest area (and vice versa), while China's afforestation projects were largely implemented in drier and previously unmanaged and unforested lands, increasing the total forest area without replacing croplands (Yu et al., 2022). Similar to net f_LULUCF, the highest emission estimates are generally derived from BLUE22 and the lowest emissions predominantly from H&N22. As stated above, this can mainly be explained by different process representation and parametrization in the models (compare Sect. 4.3). The exception of higher OSCAR22 estimates in Brazil in recent decades likely results from higher deforestation rates since 2004 and shorter turnover times for HWPs in OSCAR compared to the other BKs. In addition, OSCAR uses the averaged biome-specific carbon densities taken uniformly over the country which may overestimate emissions in particular in large countries covering differing types of the same biome (e.g. different types of forest), if land-use transitions predominantly happen in regions with lower carbon densities.

The highest differences in gross removal estimates among the BKs are found in India, Russia and the United States. In India, this may result from greater removals due to the inclusion of sub-grid-scale transitions in BLUE22 and OSCAR22, while H&N22 estimates rather negligible removals. It is noteworthy that in India, the large uncertainties in gross emissions and removals from BKs translate into a huge uncertainty in net f_LULUCF in the 1950s, but subsequently uncertainties in gross fluxes cancel out, yielding only small uncertainties for the net flux. In Russia, the models agree in a decreasing removal trend despite a considerable spread, with large removal estimates by H&N22 and small estimates by OSCAR22. In recent decades, these decreasing removals in Russia can partly be explained by the decreasing trend in the abandonment sink as was shown for the BLUE model by Winkler et al. (2023) in addition to intensified logging and wood harvest activities that cause ongoing deforestation (Kuzminyh et al., 2020). In the United States, the large BK range for net fluxes is predominantly due to large uncertainties in the removal estimates, while the gross emission estimates agree well among BKs. This removal-driven net f_LULUCF uncertainty can be explained by the inclusion of fire management in the United States in H&N22, leading to large removal estimates, while BLUE22 and OSCAR22 show much lower removal estimates.

4.5 Ratio of net-to-gross f_LULUCF from bookkeeping models

To further investigate the importance of gross fluxes, we calculate the ratio of net f_LULUCF to the sum of gross f_LULUCF (net-to-gross ratio, with the sum of gross fluxes defined as the range between gross emissions and gross removals) for the three BKs (Fig. 6; refer to Fig. A8 for RECCAP2 regions and Figs. A17 and A18 for all 186 investigated countries). Ratios close to 1 (close to −1) indicate that the net flux reflects mostly gross emissions (gross removals) and only very small gross removals (gross emissions). Ratios between 0 and 0.5 (between 0 and −0.5) indicate that the net flux represents only a small fraction of the occurring gross sources (sinks), which is the case in most countries during most of the time investigated. Ratios close to zero indicate that gross fluxes are largely compensating each other, which might indicate a sustainable land management that causes gross removals to largely offset gross emissions.

Figure 6 Decadal mean ratio of net-to-gross fluxes from land use, land-use change and forestry (with “gross” defined as the range between gross emissions and gross removals) from 1950 to 2020 derived by three bookkeeping models (BKs). The eight countries with the highest cumulative emissions since 1950 are shown in decreasing order and the United States. Values close to 1 or −1 indicate that either gross removals or gross emissions are close to 0 and the net value corresponds to either gross emissions or removals, with few compensating effects. Values near zero indicate that emissions and removals largely compensate each other. A negative ratio (green background) indicates net removals; that is, gross removals are greater than gross emissions.

Seven out of the nine investigated countries (Brazil, China, India, Canada, Russia, Nigeria and the United States) show decreasing country-level ratios of net-to-gross f_LULUCF from 1950 onward, indicating increasingly compensating gross emissions and gross removals. This is mostly due to increasing gross removals from LULUCF, while at the same time the gross emissions did not increase in such a pronounced way (particularly in Brazil, Canada, China, India and Nigeria). In some of the countries, the ratios became even negative over time (China, India, Russia, and the United States), indicating that gross removals became larger than gross emissions. Large negative net-to-gross ratios indicate that gross removals are much larger than gross emissions, and thus the (negative) net flux is mostly controlled by carbon removals from the atmosphere (e.g. in the most recent decade in European countries, Japan and Türkiye; compare Fig. 7a and Tables A1–A3).

Figure 7 Maps of (a) country-level ratio of net-to-gross carbon fluxes from land use, land-use change and forestry (with “gross” defined as the range between gross emissions and gross removals), (b) net f_LULUCF, (c) gross emissions, and (d) gross removals as average values of the three bookkeeping models for the period 2011–2021. Hatching in (a) indicates countries with a very low range in gross fluxes (average gross emissions minus gross removals smaller than 0.1 t C ha⁻¹ yr⁻¹). Green (brown) colours in (a) depict negative (positive) net fluxes with a value of $-\mathrm{1}(+\mathrm{1})$, indicating that no carbon emissions (removals) occur. The maps depict the median from three bookkeeping models (BLUE22, H&N22, OSCAR22). Grid cells with a gross flux range smaller than 0.02 tC ha⁻¹ yr⁻¹ are excluded. Gross fluxes from OSCAR22 and H&N22 were distributed using the spatial patterns of the gross flux density in BLUE for each country, respectively.

In contrast, increasing net-to-gross ratios over time are found in Indonesia and DR Congo (particularly, in the most recent decade), mainly due to strongly increasing gross emissions, which are not compensated by equally large gross removals, despite increasing removals also observable in these countries (see Fig. 5). High positive net-to-gross ratios in the most recent decade reveal large gross emissions that are not compensated by gross removals and are mainly found in the tropics and the Southern Hemisphere, in particular in Argentina, Angola, Paraguay, Bolivia, Papua New Guinea and Tanzania (Fig. 7a).

Large uncertainties in the net-to-gross ratio in Canada and China in the second half of 21st century are caused by strongly varying emission estimates and in the same period in the United States, from strongly varying removal estimates.

Near-zero net-to-gross ratios indicate that gross fluxes counterbalance each other; i.e. gross removals compensate gross emissions. Country-level near-zero net-to-gross ratios can be found in China, India, Russia and the United States, particularly in the most recent decade (compare Figs. 6 and 7a). Grid-cell-wise analysis revealed, however, that pronounced gross fluxes occurred also in these countries. This highlights that considering only the net land-use CO₂ fluxes might miss the importance of potentially large gross fluxes. This is especially true when net fluxes are estimated on a larger scale, such as at the country level and particularly for very large countries, since here the opposing gross fluxes, which often occur spatially separately, are more likely to be offset (compare, for example, United States and China in Fig. 7c and d). This, furthermore, highlights the need for spatial explicit analysis of the net as well as gross f_LULUCF and that country commitments based on net LULUCF fluxes can still be associated with large emission fluxes. Similarly, the rather vague commitment of the Glasgow Leaders' Declaration on Forests and Land Use to halt deforestation, without stating whether this accounts for gross or net transitions, can lead to strongly varying forest flux trajectories (Gasser et al., 2022). In line with Gasser et al. (2022), we therefore argue that climate mitigation measures should focus on gross fluxes from LULUCF rather than net fluxes.

5 Data availability

The NGHGI DB and the adjusted NGHGI DB data can be found under https://doi.org/10.5281/zenodo.7650360 (Grassi et al., 2023b), and OSCAR22 output can be found under https://doi.org/10.5281/zenodo.7313498 (Gasser and Shrivastav, 2022). TRENDY simulations are available via request to s.a.sitch[at]exeter.ac.uk. A dataset covering all aggregated country data for the period 1950–2021 for all datasets used in this study can be found under https://doi.org/10.5281/zenodo.8144174 (Obermeier et al., 2023).

6 Conclusions

In this study, we have comprehensively compiled country-level data on carbon fluxes from land use, land-use change and forestry (f_LULUCF) from modelling and country-report-based approaches for 186 countries. The increasing spatial resolution of modelling approaches makes it possible to provide model-based f_LULUCF estimates at country level, which can be compared to estimates based on official country reports. The comparison of multiple approaches for estimating f_LULUCF showed a fair agreement in the majority of countries, although with large differences in some other countries. The modelling approaches (BKs and DGVMs) yield generally consistent f_LULUCF estimates for the nine investigated countries. Differences, particularly across BKs, are due to differences in land-use forcing data, process implementation and parametrization.

Similarly, DGVM estimates strongly depend on the land-use forcing data. For some of the investigated countries, further uncertainties in the DGVM estimates of a similar magnitude are caused by the environmental forcing data used, namely present-day environmental forcing, which is more comparable to BKs and country-report-based approaches, compared to transient environmental forcing, which better reflects historical environmental changes in the real world. In the majority of investigated countries, f_LULUCF estimates based on official country reports (NGHGI DB and FAOSTAT) are lower compared to the modelled estimates. However, once the varying characteristics and definitions (in particular the so-called managed land proxy) are accounted for, the differences become substantially lower in most countries.

Analysing the gross fluxes from BKs revealed that short-term variations in net fluxes are mostly linked to gross emissions, which show large temporal variability, while gross removals rather impact the long-term trends of net fluxes. Uncertainties in net fluxes mainly relate to uncertainties in gross emissions (for example, in Brazil, Canada, China and DR Congo) but can also be strongly impacted by uncertainties in gross removals (for example in the United States). In India and Russia, pronounced uncertainties in both gross emissions and gross removals largely compensate each other, which results in rather low uncertainties of the net flux. Furthermore, the investigation of the net-to-gross ratio revealed that the net flux is comprised by large gross fluxes in most countries and over most of the investigated time from 1950 onward. It is noteworthy that gross fluxes increasingly compensate each other in most of the countries over time. Considering only net fluxes might thus miss potentially important and large gross fluxes. In addition, grid-cell-wise analysis revealed that pronounced gross and net fluxes may occur within a country at different locations, even though net fluxes are close to zero when averaged at the country level, which highlights the need for spatially explicit data on gross fluxes.

Consequently, model-based spatial data as presented in this study may support the identification of component-wise, historical and/or regional “uncertainty hotspots” that particularly need improved f_LULUCF estimations. For example, the uncertain emission estimates from BKs in India, from the 1980s onward in China, and for the most recent decade in Brazil and DR Congo could be improved by better land-use forcing data. Likewise, the uncertain removal estimates in Russia and India could be improved by incorporating better land-use forcing data and improved process representation in models, such as for fluxes related to land abandonment. The differences in f_LULUCF estimates in these “hotspots” highlight the need for a careful interpretation of the outputs from the varying methods and for a further reconciliation of the different approaches, in particular regarding the different components considered and the methods used for their estimation. For example, in Canada and Nigeria, increasing differences between modelled estimates and those from the NGHGI DB, even after adjustment of the latter, call for in-depth analysis of the underlying drivers.

To this end, we argue for a systematic model evaluation and improved parametrization of models, in particular regarding land-use forcing data, parametrized carbon densities, and the different processes represented in the models. In addition, the definition and framework issues implicitly underlying all datasets that still exist, for example, definition/inclusion of LASC, whether there are transient C densities or not, biome and plant functional type (PFT) definitions, and managed land proxy, should be addressed. To further increase the confidence in f_LULUCF estimates, more approaches similar to bookkeeping models, for example, by DGVM simulations under BK-like protocol, could be used. In addition, all models (BKs and DGVMs) should use more spatially explicit forcing datasets. In particular, Earth observation data may provide improved spatially explicit data (e.g. of land degradation and restoration efforts) on carbon densities in vegetation and soil, forest regrowth by incorporation of forest age classes, and forest management from optical and microwave satellite measurements, as well as data on carbon fluxes from atmospheric inversions. Furthermore, spatially explicit data on land-use activities from Earth observations could improve the quality of report-based estimates where data sources are scarce and improve comparability with estimates from modelling approaches. Such an improved incorporation of Earth observation data into modelling and country-report-based approaches may provide substantial advancements in the assessment and understanding of CO₂ fluxes from LULUCF.

Appendix A: Description of individual approaches for estimating f_LULUCF

A1 Bookkeeping models

BKs use spatial information on land-use activities to derive net f_LULUCF by summing up all gross carbon fluxes that occur due to land conversion and land management (Pongratz et al., 2014). To estimate carbon fluxes from LULUCF, BKs rely on observation-based carbon stock densities, growth curves (uptake) and decomposition curves (release) of soil and vegetation carbon that are specific for each conversion type. Fluxes due to land-use conversion can occur instantaneously (fluxes upon or within the year of LULUCF) and in the years following the conversion (legacy fluxes, for example, from readjustment of carbon stocks).

Figure A1 Regions as defined by REgional Carbon Cycle Assessment and Processes Phase 2 (RECCAP2; Tian et al., 2019).

Figure A2 Net per-capita carbon fluxes from land use, land-use change and forestry (f_LULUCF) from three bookkeeping models (BKs; data from GCB2022 simulations). (a) Cumulative per-capita carbon fluxes over 1950–2021 and (b) average per-capita carbon fluxes in 2011–2021. The bars show the mean of the three BKs (filled bars) and minimum and maximum estimates (hatched bars). Numbers in parentheses show the multi-model average and standard deviation (in tC per capita in (a) and tC per capita yr⁻¹ in (b)). Colours indicate the absolute quantities, showing countries with net emissions in red and countries with net removals in green. All 186 country aggregates from this study are shown in decreasing order of their (a) cumulative and (b) most recent annual f_LULUCF. In each panel, the top 10 emitters and the five countries with the largest removals are labelled. The figure corresponds to Fig. 1 in the main paper, with the difference that fluxes are shown per capita. Note that values for very small countries should be interpreted with care as the relatively low resolution of many models creates uncertainty at the small scale.

Figure A3 Net per-area carbon fluxes from land use, land-use change and forestry (f_LULUCF) from three bookkeeping models (BKs; data from GCB2022 simulations). (a) Cumulative per-area carbon fluxes over 1950–2021 and (b) average per-area carbon fluxes in 2011–2021. The bars show the mean of the three BKs (filled bars) and minimum and maximum estimates (hatched bars). Numbers in parentheses show the multi-model average and standard deviation (in tC ha⁻¹ in (a) and tC ha⁻¹ yr⁻¹ in (b)). Colours indicate the absolute quantities, showing countries with net emissions in red and countries with net removals in green. All 186 country aggregates from this study are shown in decreasing order of their (a) cumulative and (b) most recent annual f_LULUCF. In each panel, the top 10 emitters and the five countries with the largest removals are labelled. The figure corresponds to Fig. 1 in the main paper, with the difference that the fluxes are shown per area. Note that values for very small countries should be interpreted with care as the relatively low resolution of many models creates uncertainty at the small scale.

Figure A4 Time series of net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2021 derived by bookkeeping models (BKs) and TRENDYv8 simulations with dynamic global vegetation models (DGVMs; used in the GCB2019) under present-day climate forcing for the RECCAP2 regions. Regions are sorted according to their cumulative net emissions from 1950–2021, as derived by three bookkeeping models. The complete region designations for the acronyms used can be found in the legend of Fig. A1. The figure shows the mean and absolute range of three BKs (using the GCB2022 simulations, BKs 2022) and the median and interquartile range of the eight DGVMs. Additionally, estimates from BLUE simulations from the GCB2019 (BLUE19; blue solid) and from the GCB2022 (BLUE22; blue dashed) are shown to illustrate the impact of updates in the LUHv2 forcing data. BLUE19 data are only available until 2019 and TRENDYv8 data are only available until 2018. Greenish background depicts negative f_LULUCF, which is carbon removal from the atmosphere. The figure corresponds to Fig. 2 in the main paper, with the difference that the fluxes are shown for the RECCAP regions.

Figure A5 Time series of net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2021 derived from TRENDYv8 simulations with dynamic global vegetation models (DGVMs) under historical (transient) and fixed present-day environmental conditions (compare Appendix) for the RECCAP2 regions. Regions are sorted according to their cumulative net emissions from 1950–2021, as derived by three bookkeeping models. The complete region designations for the acronyms used can be found in the legend of Fig. A1. The figure shows the median and interquartile range of the eight DGVMs. Greenish background depicts negative f_LULUCF, which is carbon removal from the atmosphere. The figure corresponds to Fig. 3 in the main paper, with the difference that the fluxes are shown for the RECCAP regions.

Fluxes from peat fire and peat drainage are not directly modelled by BKs but are added from external data. For the GCB2022 simulations, which we use here, peat fire emissions were added from the Global Fire Emission Database (GFED4s; van der Werf et al., 2017) for Brunei Darussalam, Indonesia, Malaysia and Papua New Guinea. Peat drainage emissions were added for all countries as the average from FAO data (Conchedda and Tubiello, 2020) and DGVM simulations with ORCHIDEE-PEAT (Qiu et al., 2021) and LPX-Bern (Lienert and Joos, 2018; Müller and Joos, 2021). More details can be found in Friedlingstein et al. (2022a).

Figure A6 Net carbon flux from land use, land-use change and forestry (f_LULUCF) in 1950–2019 derived from individual bookkeeping models (BKs) used in GCB2022 (BKs 2022) and country-report-based estimates (FAOSTAT and NGHGI DB) for the RECCAP2 regions. The use of dashed and solid lines indicates that BK estimates and country-report-based estimates are not directly comparable (see Sect. 2.5 and the Appendix). For better comparability, the adjusted NGHGI DB estimates, matching the f_LULUCF definition of the BKs, are additionally shown. The grey line (shading) depicts the median (range) of the three BKs. The light-green background indicates negative f_LULUCF that is net carbon removal from the atmosphere by LULUCF. The figure corresponds to Fig. 4 in the main paper, with the difference that the fluxes are shown for the RECCAP regions.

The three BKs notably differ in (1) the implemented processes; (2) the land-use forcing data; and (3) the assigned carbon densities, response curves and pool allocation fractions (compare Friedlingstein et al., 2022b; Bastos et al., 2021a):

1. BLUE22 and OSCAR22 include sub-grid-scale transitions between all vegetation types (e.g. from shifting cultivation – a rotation cycle between forest and agriculture), which presumably leads to higher emissions, whereas sub-grid-scale transitions are not implemented in H&N22 – only if a country's forest loss reported to FRA exceeds agricultural expansion based on FAO land-use data does H&N22 assume that this area is cleared for shifting cultivation. BLUE22 includes gross fluxes related to degradation from primary to secondary land in the case that natural vegetation is used as rangeland. H&N22 considers fire management in the United States and southeast Asia, in contrast to BLUE22 and OSCAR22.

2. BLUE22 uses the LUH2-GCB2022 dataset (an update to most recent harmonized land-use change data (LUH2 v2h); Chini et al., 2021; Hurtt et al., 2020) based on HYDE3.3, whose contemporary land use is constrained by annual ESA CCI Land Cover and updated agricultural areas from FAO (Klein Goldewijk et al., 2017). These data are globally consistent but have a relatively coarse spatial resolution (0.25^∘ × 0.25^∘) and may thus exclude regional and local specifics (Bastos et al., 2018; Li et al., 2018; Kondo et al., 2022). H&N22 uses FAO-FRA data from 2020 for forest (from 1990 onward and various sources before; FAO, 2020) and data from FAOSTAT for other land uses and applies a 5-year running mean on the activity data before flux calculations (Friedlingstein et al., 2022b). OSCAR22 calculates a best-guess estimate of f_LULUCF based on a combination of the LUH2-GCB2022 dataset and FAO-FRA data (Gasser et al., 2020, 2022). BLUE22 output is spatially explicit, while H&N22 and OSCAR22 provide country-level estimates.

   We further use BLUE data from GCB2019 (hereafter BLUE19), forced with LUH2-GCB2019 data based on HYDE3.2 (Klein Goldewijk et al., 2011; Chini et al., 2021). As the BLUE model code was not changed between GCB2019 and GCB2022, this allows the impact of changes in the LULUCF forcing data to be isolated, with HYDE3.2 using 1 year of ESA CCI as a reference year for the spatial land cover patterns, whereas HYDE3.3 uses time-varying ESA maps, which led to spatio-temporally improved land cover maps (Rosan et al., 2021). From 2018 onward, the underlying LUH2 data linearly interpolate the trend in cropland, pasture and urban area of the previous 5 years until the year 2021. This approach does not properly reflect the dynamics in regions with intensive land use and land-use changes (LULUCs) in most recent years; therefore, the LUH2-GCB2022 land-use forcing data in Brazil are taken from the MapBiomas dataset (collection 6) for the period 1985–2020, according to the approach described in Friedlingstein et al. (2022a).

3. H&N22 assigns vegetation carbon densities at country level (based on official country reports) and soil carbon densities globally for 20 types of ecosystems (Houghton and Nassikas, 2017). BLUE22 assigns vegetation and soil carbon densities for 11 PFTs globally (based on literature values; Houghton et al., 1983; Hansis et al., 2015). OSCAR22 uses carbon densities derived from DGVMs and response curves specified for 96 world regions and five biomes (Gasser et al., 2020, 2022). Similarly, carbon response curves in H&N22 and BLUE22 are assigned for each of the 20 ecosystems and 11 PFTs considered, respectively. The BLUE model, in general, has highest carbon densities implemented which causes high f_LULUCF estimates, as was shown by Bastos et al. (2021a, b), where parameterizing BLUE with H&N carbon densities led to a 24 % reduction of global cumulative f_LULUCF from 1850–2015. Using carbon densities and response curves that are static over time, BLUE22 and H&N22 do not explicitly model the effects of environmental changes (e.g. increased CO₂ concentration and climatic change), although some are implicitly captured within the observed carbon densities and response curves (Pongratz et al., 2014). In contrast, OSCAR22 includes transient environmental response due to its calibration to transient DGVM simulations (Gasser et al., 2017, 2020). In addition, the allocation of HWPs to different product pools differs between the models: BLUE22 uses three HWP pools (with turnover times of 1, 10, and 100 years) with fixed allocation fractions for each PFT. H&N22 assigns time-variant fractions for five pools (fuel and industrial and 1-year, 10-year and 100-year turnover times) specific for each country (Bastos et al., 2021a). OSCAR22 uses three HWP pools (with average turnover times of 0.75, 6.0 and 65 years) and allocation fractions specific to regions and biomes (Gasser et al., 2017, 2020).

In some countries depicted in the Appendix, particularly in arid world regions, gross emissions from BKs are negative in some years. This relates to the definition of gross emissions which include emissions from deforestation, forest degradation and wood harvest and fluxes from transitions between natural land, cropland and pasture. The latter, however, may cause negative carbon fluxes if carbon densities of the initial land cover are lower than the carbon densities of the converted land cover. Consequently, this may cause gross emissions to be negative, particularly in dryland countries with little forest cover. On the global scale, this effect is negligible.

Figure A7 Gross and net fluxes from land use, land-use change and forestry (f_LULUCF) from 1950 until 2021, as derived by three bookkeeping models (BKs; used in the GCB2022) for the RECCAP2 regions. Regions are sorted according to their cumulative net emissions from 1950–2021, as derived by three bookkeeping models. The complete region designations for the acronyms used can be found in the legend of Fig. A1. Greenish background depicts negative f_LULUCF, which is carbon removal from the atmosphere. The figure corresponds to Fig. 5 in the main paper, with the difference that the fluxes are shown for the RECCAP regions.

Figure A8 Gross and net fluxes from land use, land-use change and forestry (f_LULUCF) from 1950 until 2021, as derived by three bookkeeping models (BKs; used in the GCB2022) for the RECCAP2 regions. Regions are sorted according to their cumulative net emissions from 1950–2021, as derived by three bookkeeping models. The complete region designations for the acronyms used can be found in the legend of Fig. A1. Lines (shaded area) depict the mean (range) of the BK estimates. The figure corresponds to Fig. 6 in the main paper, with the difference that the fluxes are shown for the RECCAP regions.

Figure A9 As in Fig. 2 but for the 93 top-emitting countries according to their cumulative emission in 1950–2021, as derived by three bookkeeping models (without the lines for the differences of the estimations; remaining countries are shown in Fig. A10). Green-shaded areas depict the range of net carbon removals. For the complete country designations of the alpha-3 codes, refer to Tables A1–A3.

A2 Dynamic global vegetation models

DGVMs are used in the GCB for the uncertainty assessment of f_LULUCF and the estimation of the natural land sink. Moreover, DGVMs are frequently used in more detailed studies, for example, on the effects of land-use changes on local climate, due to their implementation of complex biogeochemical and biogeophysical processes and their capacity to simulate transient environmental responses (e.g. Krause et al., 2018; Winckler et al., 2017; Bright et al., 2017). We aggregated country-level f_LULUCF based on DGVMs that performed simulations within the project “Trends and drivers of the regional-scale emissions and removals of carbon dioxide” (TRENDY; Le Quéré et al., 2014; Sitch et al., 2015).

Figure A10 As in Fig. A9 but showing the lowest-emitting countries and countries with net removals according to their cumulative fluxes in 1950–2021, as derived by three bookkeeping models. Green-shaded areas depict the range of net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A11 As in Fig. 3 but for the 93 top-emitting countries according to their cumulative emission in 1950–2021, as derived by three bookkeeping models (without the lines for the differences of the estimations; remaining countries are shown in Fig. A12). Green-shaded areas depict the range of net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

DGVMs do not directly output f_LULUCF. Instead, f_LULUCF is estimated as the difference between two simulations with the same environmental forcing, one including and one excluding LULUCF. To derive f_LULUCF for each grid cell and each (yearly) time step, the net biome productivity (NBP) of the simulation including LULUCF is subtracted from the NBP in the simulation excluding LULUCF, the latter using a constant pre-industrial LULUCF map (from 1700) over time. Thereby, f_LULUCF from DGVMs includes instantaneous as well as legacy fluxes, like BK estimates.

Figure A12 As in Fig. A11 but showing the lowest- emitting countries and countries with net removals according to their cumulative fluxes in 1950–2020, as derived by three bookkeeping models. Green-shaded areas depict the range of net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A13 As in Fig. 4 but for the 93 top-emitting countries according to their cumulative emission in 1950–2020, as derived by three bookkeeping models (remaining countries are shown in Fig. A14). Green-shaded areas depict the range of net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A14 As in Fig. A13 but showing the lowest-emitting countries and countries with net removals according to their cumulative fluxes in 1950–2020, as derived by three bookkeeping models. Green-shaded areas highlight net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A15 As in Fig. 5 but for the 93 top-emitting countries according to their cumulative emission in 1950–2021, as derived by three bookkeeping models (remaining countries are shown in Fig. A15). Green-shaded areas highlight net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A16 As in Fig. A15 but showing the lowest-emitting countries and countries with net removals according to their cumulative fluxes in 1950–2021, as derived by three bookkeeping models. Green-shaded areas highlight net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

DGVM simulations can be forced with different environmental conditions, where some environmental variables are set constant (fixed) at either pre-industrial or present-day levels or follow observed, transient conditions (for an overview of the different simulations, refer to Obermeier et al., 2021). Here, we use simulations under fixed present-day environmental conditions as they most closely resemble BK simulations and country-report-based approaches (using observed C densities) and are recommended by RECCAP2 (Ciais et al., 2022). Present-day environmental simulations are run under the CO₂ concentration from 2018 throughout the simulated period and recycle the climate from 1999–2018 using the mean and variability of the individual years in this period (Obermeier et al., 2021). DGVM simulations with present-day conditions are not performed every year. The most recent present-day simulations available stem from TRENDYv8 used in GCB2019 (Friedlingstein et al., 2019).

We additionally employ transient DGVM simulations from TRENDY v8 as transient simulations are operationally available, commonly used within the scientific community, and enable us to derive the difference in f_LULUCF that results from different environmental forcing (by comparing to present-day TRENDYv8 simulations). The transient environmental simulations used here are forced with observation-based temperature, precipitation and incoming surface radiation data (0.5^∘ × 0.5^∘ resolution) of the Climatic Research Unit (CRU) and the Japanese Reanalysis (JRA; Friedlingstein et al., 2019; Harris et al., 2014) and global atmospheric CO₂ concentrations from ice core data before 1958 (Joos and Spahni, 2008) combined with National Oceanic and Atmospheric Administration (NOAA) data from 1958 onward (Dlugokencky and Tans, 2020).

To enable a robust comparison between the different forcings, we selected only those nine DGVMs that provide present-day in addition to transient simulations within TRENDYv8: CLASS-CTEM (Melton and Arora, 2016), DLEM (Tian et al., 2015), JSBACH (Mauritsen et al., 2019), LPJ-GUESS (Smith et al., 2014), LPX-Bern (Lienert and Joos, 2018), ORCHIDEE (Krinner et al., 2005), ORCHIDEE-CNP (Goll et al., 2017), SDGVM (Walker et al., 2017) and VISIT (Kato et al., 2013).

Differences among the DGVMs mainly result from (1) differing model parameters and (2) the implementation of different processes with varying complexities (Friedlingstein et al., 2022b):

1. Model parametrization differs, for instance, in the distinction of primary and secondary forests and turnover rates of product pools (Kondo et al., 2022), the fraction of directly emitted carbon upon LULUCs, and the implemented decomposition rates and resulting soil carbon densities (Goll et al., 2015).

2. Implemented processes differ, for example, as some DGVMs include fires (without distinguishing whether they are natural or anthropogenic), while others have no fire implemented. Other natural disturbances are not included by the DGVMs. Additionally, some DGVMs consider nitrogen or phosphorus cycles, cropland irrigation, shifting cultivation, forest degradation, and residue carbon after deforestation and wood and crop harvests, while others do not (for more details refer to Bastos et al., 2020b; Kondo et al., 2022; Friedlingstein et al., 2022b). It is noteworthy that the inclusion of processes such as shifting cultivation, wood harvest, grazing, crop harvest and cropland management increases historic (1901–2014) LULUCF emissions from DGVMs by 20 %–30 % for each of these processes (Arneth et al., 2017).

A3 National greenhouse gas inventories under the UNFCCC

To take stock and track progress towards their NDCs, countries report official inventory statistics of GHG emissions and removals to the UNFCCC via different schemes. In line with Grassi et al. (2023a), in this study, we refer to any of such official country reports on anthropogenic GHG data submitted to UNFCCC as national greenhouse gas inventories (NGHGIs).

Following the reporting guidelines of UNFCCC, NGHGIs are submitted regularly (annually for Annex I countries and typically at least every few years for non-Annex I countries). NGHGIs are required to meet the key principles of transparency, accuracy, completeness, consistency and comparability (TACCC). Nonetheless, depending on national circumstances and conditions, reporting is done with different frequency and sophistication of the underlying methods, and large methodological uncertainties exist (Federici et al., 2017). According to the IPCC best-practice guidelines (IPCC, 2006), some flexibility is allowed for forest definitions within the NGHGIs; for example, thresholds of parameters defining forests (minimum area, tree crown cover, tree height) can be chosen from an allowed interval (0.05–1.0 ha², 10 %–30 %, 2–5 m, respectively) in the first national communication but must be kept constant afterwards. The IPCC defined different tiers that indicate the sophistication of the applied methods ranging from Tier 1 to Tier 3, with Tier 3 being the most demanding in terms of data availability and method complexity. Most Annex I countries (members of the OECD and some transitional economies) report all land-use fluxes annually since 1990 following the 2006 IPCC Guidelines (and partly following the 2019 IPCC refinements) (IPCC, 2006, 2019). Guidelines for non-Annex I countries (emerging economies, for example, Brazil, China, Democratic Republic of the Congo, India, Indonesia and Nigeria) are more flexible, the applied methodologies are generally less complex and reporting often started in 2000 only (IPCC, 2006).

Table A1 Mean, range, minimum and maximum estimates of country-level annual mean (2011–2021) and cumulative (1950–2021) net f_LULUCF, as derived from three bookkeeping models, and the respective country's rank. Countries are sorted in alphabetical order. Note that the rank of the cumulative f_LULUCF (in bold) gives the order of countries as they occur in Figs. A9–A18.

Table A2 Mean, range, minimum and maximum estimates of country-level annual mean (2011–2021) and cumulative (1950–2021) net f_LULUCF, as derived from three bookkeeping models, and the respective country's rank. Countries are sorted in alphabetical order. Note that the rank of the cumulative f_LULUCF (in bold) gives the order of countries as they occur in Figs. A9–A18.

Table A3 Mean, range, minimum and maximum estimates of country-level annual mean (2011–2021) and cumulative (1950–2021) net f_LULUCF, as derived from three bookkeeping models, and the respective country's rank. Countries are sorted in alphabetical order. Note that the rank of the cumulative f_LULUCF (in bold) gives the order of countries as they occur in Figs. A9–A18.

The most important developing countries (Brazil, Indonesia, China, India and Mexico) rely on national inventories and, in the case of Brazil, even have a NGHGI comparable to some of the developed countries. In contrast, NGHGIs from most other non-Annex I countries rely on empirical emission factors to estimate f_LULUCF, which are representative rates of emissions (e.g. for specific forest and climate types) that are usually obtained from averaged measurement data sampled under certain environmental conditions and, thus, hardly capture local dynamics. Such basic approaches using default values correspond to Tier 1 level, while reports from Annex I countries are often based on country-specific statistical or process-based models using national activity data (i.e. Tier 2 or Tier 3 level).

To reduce the uncertainties associated with the different reporting schemes, we use the newly compiled data from Grassi et al. (2023a), which is based on official country reports to the UNFCCC but additionally includes quality checks and gap-filling if necessary (hereafter named NGHGI DB). NGHGI DB data are freely available under https://doi.org/10.5281/zenodo.7650360 (Grassi et al., 2023b) from 2000 onward and are considered a realistic approximation of the data to be used for the forthcoming Global Stocktake in 2023 (Grassi et al., 2023a).

NGHGIs are supposed to encompass all LULUCF fluxes from areas considered managed, including forest land, cropland, grassland, wetlands and settlements, as well as emissions from organic soils and fires (and some also include shifting cultivation). Thereby, the reports should include pools for dead wood, litter, soil organic C and HWPs. However, many non-Annex I countries report only fluxes from deforestation, and only few include fluxes from other LULUCF categories. To distinguish between anthropogenic and non-anthropogenic fluxes, national authorities use the “managed land” proxy where emissions and removals on managed lands are counted for, whereas fluxes on unmanaged lands are not reported.

A4 FAOSTAT

The FAOSTAT f_LULUCF data are a component of the FAO Emissions database, developed to assess the role of food and agriculture in global anthropogenic emissions (Tubiello, 2019; Tubiello et al., 2022; FAO, 2021). The FAOSTAT f_LULUCF data cover carbon emissions and carbon removals as well as non-CO₂ emissions from biomass fires (not considered in this work) on the following IPCC land-use and land-use-change categories: (1) emissions and removals on forests and from deforestation (Tubiello et al., 2021), (2) emissions from peatland fires (Rossi et al., 2016; Prosperi et al., 2020), and (3) emissions from peatland drainage (Conchedda and Tubiello, 2020; Tubiello et al., 2016).

FAOSTAT f_LULUCF data cover 238 countries and territories, with sub-regional, regional and global aggregates for the period 1990–2020. They are estimated by applying IPCC (2006) guidelines to LULUCF activity data generated either through official country reporting processes, such as the FRA (FAO, 2020), or through analysis of geospatial data carried out by FAO under its mandate. Specifically, forest fluxes are based on carbon stock change statistics computed directly at national level, following carbon stock and (net) area change statistics reported by member countries to FAO at 5-year intervals (FAO, 2020; Tubiello et al., 2021). Deforestation emissions are computed separately for the two FAO forest types, naturally regenerating forest (comprising primary and secondary natural forests) and planted forest. It is noteworthy that FAOSTAT estimates only include the above- and below-ground biomass pools, while fluxes resulting from changes in other pools, such as the soil carbon pools, are not modelled, except for those in organic soils. The estimated emissions and removals on forest land do not distinguish between anthropogenic and non-anthropogenic fluxes and include indirect climate and CO₂ effects in line with the IPCC management land proxy. Additionally, FAOSTAT f_LULUCF data comprise anthropogenic peat fires computed at grid cell level using the map of Histosols from the Harmonized World Soils Database (as a proxy for spatial peatland distribution) combined with remote sensing products of MODIS burnt area and underlying MODIS land cover maps. IPCC (2013) tier 1 emission factors are associated with each land cover/ecological zone by merging with the IPCcarbon-JRC agro-climatic zone map (FAOSTAT, 2020), and results are aggregated at country level (for details, refer to Rossi et al., 2016; Prosperi et al., 2020). Similarly to the BKs, peat fire emissions are only considered anthropogenic for Southeast Asian countries, albeit data are available globally. Emissions from drained peatlands, also computed at grid cell level, are estimated based on the harmonized world map of Histosols and the ESA CCI land cover map to identify cropland areas, assuming cultivation on peatland area is a proxy for anthropogenic drainage (Conchedda and Tubiello, 2020).

Figure A17 As in Fig. 6 but for the 93 top-emitting countries according to their cumulative emission in 1950–2021, as derived by three bookkeeping models (remaining countries are shown in Fig. A18). Green-shaded areas highlight net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.

Figure A18 As in Fig. A17 but showing the lowest-emitting countries and countries with net removals according to their cumulative fluxes in 1950–2020, as derived by three bookkeeping models. Green-shaded areas highlight of net carbon removals. For the official country names and each country's rank (order of subplots), refer to Tables A1–A3.