Uncertainties in carbon emissions from land use and land cover change in Indonesia

Indonesia is currently one of the three largest contributors of carbon emissions from land use and land cover change (LULCC) globally, together with Brazil and the Democratic Republic of the Congo. However, until recently, there were only limited reliable data available on LULCC across Indonesia, leading to a lack of agreement on the drivers, magnitudes, and trends of carbon emissions between different estimates. Accurate LULCC should improve robustness and reduce the uncertainties of carbon dioxide (CO₂) emissions from land use change (ELUC) estimation. Here, we assess several cropland datasets that are used to estimate ELUC in dynamic global vegetation models (DGVMs) and bookkeeping models (BKMs). Available cropland datasets are generally categorized as either statistics-based, such as the Food and Agricultural Organization (FAO) annual statistical dataset, or satellite-based, such as the Mapbiomas dataset, which is derived from Landsat satellite images. Our results show that national-statistics-based and satellite-based estimates have little agreement on temporal variability and cropland area changes. On some islands, they show spatial similarity, but differences appear on the main islands such as Kalimantan, Sumatra, and Java. These differences lead to spatiotemporal uncertainty in carbon emissions. The different land cover forcings (national-statistics-based and satellite-based) in a single model (JULES-ES) result in ELUC uncertainties of about 0.08 [0.06 to 0.11] PgC yr⁻¹. Furthermore, we found that uncertainties in ELUC estimates are also due to differences in the carbon cycle models in DGVMs, as DGVMs driven by the same land cover dataset show differences in ELUC estimates of 0.12 ± 0.02 PgC yr⁻¹ with a 95 % confidence level and range [−0.04 to 0.35] PgC yr⁻¹. This is consistent with other products such as BKMs that estimate 0.14 [0.12 to 0.15] PgC yr⁻¹, with both having steady trends. We also compare the emissions with those from the National Greenhouse Gas Inventory (NGHGI) product. The NGHGI estimates (based on BUR3, the periodic official government report on greenhouses gases to UNFCCC) have much lower carbon emissions (0.06 ± 0.06 PgC yr⁻¹), though with an increasing trend. These numbers double when we include emissions from peat fire and peat drainage: the DGVM ensemble indicates emissions of 0.23 ± 0.05 PgC yr⁻¹, and BKMs indicate emissions of 0.24 ± 0.01 PgC yr⁻¹. In contrast, emissions based on the Indonesian NGHGI remain much lower (BUR2: 0.18 ± 0.07 PgC yr⁻¹; BUR3: 0.13 ± 0.10 PgC yr⁻¹). Furthermore, emission peaks occur in years of moderate to strong El Niño events. Several improvements might reduce uncertainties in carbon emissions from LULCC in Indonesia, such as a combination of a satellite-based dataset with a national-statistics-based dataset, inclusion of peat-related emissions in DGVMs, and potentially explicit inclusion of palm oil in models, as this is a major crop in Indonesia. Overall, the analysis shows that carbon emissions have no decreasing trend in Indonesia. Therefore, deforestation and forest fire prevention remain vital for Indonesia.

1 Introduction

Indonesia plays a significant role in the global carbon cycle through carbon emissions from land use and land cover change (LULCC) (Friedlingstein et al., 2023). Together with Brazil and the Democratic Republic of Congo, these three countries are the highest carbon emitters in the world from land use change. Their emissions combined are more than half of the global total land use emissions (0.75 [0.59–0.87] PgC yr⁻¹ or 56 [44–68] %) (Friedlingstein et al., 2023).

Indonesian forests are spread across several islands. The major islands are Sumatra, Kalimantan, Sulawesi, and Papua (Gaveau et al., 2021). In recent decades, two islands, Sumatra and Kalimantan, have experienced extensive land use changes (Gaveau et al., 2014, 2016; Margono et al., 2012). There were many forested areas that have now been converted into agricultural lands, predominantly palm oil plantations. Forest area declined by 9.79 million hectares (Mha) or 11 % between 2001 and 2019, with near-equal expansion of oil palm plantations of 8.48 Mha (Gaveau et al., 2022). This forest loss is also confirmed by the Food and Agricultural Organization (FAO), which shows the decline from 101.12 to 92.13 Mha for the last 2 decades, or 9.14 Mha lost (FAO, 2022).

These high rates of land conversion in regions of carbon-dense and often pristine natural forests result in high carbon emissions (Hong et al., 2021). Furthermore, there are large peat and forest fires during extreme dry seasons in Indonesia (Fernandes et al., 2017; Nechita-Banda et al., 2018). This has added to the carbon released into the atmosphere from Indonesia, together with fires induced by humans (Brasika, 2023; FAO, 2023; van der Werf et al., 2017). Although there is a declining trend in reported forest fires between 2003 and 2018 (van Wees et al., 2021), fires continue to contribute significant emissions. For example, in 2019–2020, deforestation fires in Indonesia continued contributing to high greenhouse gas (GHG) emissions, i.e. about 3.7 ± 0.4 Gt CO₂ eq., and half of this can be attributed to emissions from peatland fires (Datta and Krishnamoorti, 2022).

The emission of carbon from LULCC caused by deliberate human activities is represented as emissions from land use change (ELUC) (PgC yr⁻¹) in the Global Carbon Budget (GCB) calculated using a set of bookkeeping models (Friedlingstein et al., 2023). This includes CO₂ fluxes from deforestation, afforestation, logging, forest degradation, shifting cultivation, and regrowth of forests. ELUC can also be calculated using dynamic global vegetation models (DGVMs), which are part of the Trends in Net Carbon Exchange (TRENDY) project (Friedlingstein et al., 2019; Sitch et al., 2015), which supplies ensemble DGVM results each year to the GCB. The LULCC input for the GCB is from the Land Use Harmonization (LUH2) dataset (Chini et al., 2021; Hurtt et al., 2020), which has temporal and spatial challenges (Rosan et al., 2021) when employed for national-scale analysis.

With the development of multi-year land cover datasets from Earth observations, there is potential to reduce the uncertainty in carbon emissions from LULCC. This was previously applied in the GCB 2022 (Friedlingstein et al., 2022), where the authors replaced the FAO land use data with satellite-based Mapbiomas data for Brazil (Rosan et al., 2021). This resulted in improvements in LULCC datasets that gave consistent estimates of ELUC trends using different methodologies for Brazil.

In this research, we apply a similar approach to improve our understanding of carbon emissions from LULCC in Indonesia over the last 2 decades (2000–2022). However, we focus more on the uncertainties, i.e. not only the causes by the input, but also the different approaches of models and/or inventories. We used LUH2-GCB2022, which is the harmonization of several datasets, including satellite data with FAO land use data and satellite-based cropland (Mapbiomas Indonesia 1.0 – MB1 – and Mapbiomas Indonesia 2.0 – MB2) datasets, as input to JULES-ES. First, we compare the temporal and spatial variability of cropland from these datasets. Second, we simulate ELUC with two approaches: various LULCC driver datasets (LUH2-GCB2022, MB1, and MB2) in one model (JULES-ES) and comparison with results using one LULCC driver (MB1) in various models in TRENDY v12. Last, we compare our estimates with other models or products, i.e. bookkeeping models (BKMs), the National Greenhouse Gas Inventory (NGHGI), and FAO GHG for Indonesia.

Despite addressing many uncertainties caused by the input data, model structure, and different products, we emphasize the agreement between those models and products that do appear. This includes the similarity of the trends, their variability, their magnitudes, and their potential connection to interannual variability. This might result in more robust understanding of carbon emissions from Indonesia as these agreements appear in the uncertainties. This is essential as the Indonesian contribution to global carbon emissions from land use change is massive, but this multi-dataset assessment has not been done before.

2 Methods

2.1 Land use dataset

2.1.1 Mapbiomas Indonesia dataset

Mapbiomas is an initiative aimed at mapping LULCC using a fast, reliable, and cost-effective methodology. Initially launched in Brazil in 2015, it has since been applied to other countries. Indonesia became the first country not in Latin America to adopt the method, receiving training from experts in Brazil (experts from the Asian country were trained by the Brazilian Mapbiomas team to launch the land use mapping platform; Silva et al., 2022). Mapbiomas utilizes Landsat images and employs a pixel-based approach. Machine learning algorithms have been developed to analyse these images, leveraging the powerful cloud processing capabilities provided by Google Earth Engine (GEE).

The resulting dataset, MB1, classifies the LULCC from 2000 to 2019 (Mapbiomas Indonesia, 2022). While the classification scheme is adapted from Mapbiomas Brazil, it also includes specific characteristics relevant to Indonesia. For instance, agriculture in Mapbiomas Indonesia includes a palm oil category, reflecting the significant impact of palm oil production on LULCC in Indonesia (Mapbiomas Indonesia, 2022). However, there is no specific grazing land category in Mapbiomas Indonesia, as grazing land is considerably small and steady in Indonesia. Nevertheless, we still use the grazing land from LUH2-GCB2022 for completeness (Fig. 1). In case the total combination of the Mapbiomas agriculture fraction and the LUH2-GCB2022 grazing land fraction is higher than 1 in a given grid cell, we reduce the grazing land fraction accordingly. So, the total area of cropland + grazing land will be equal to or less than 1.

Mapbiomas has a high spatial resolution of 30 m, which is significantly higher than other products such as ESA CCI LC at 300 m. The Mapbiomas dataset has the potential to detect even small changes in LULCC. This dataset finds various applications, including the calculation of CO₂ emissions resulting from land use change (Garofalo et al., 2022). Currently, Mapbiomas Indonesia is on its Collection 2.0 version (MB2). This recent dataset has completed an accuracy assessment and validation process (using ∼ 12 967 samples), with a resulting overall accuracy of 77.2 %. This MB2 covers a longer period, 2000–2022, with some improvement in land use classification. The agriculture is divided into more categories: palm oil, rice paddy, pulpwood plantation, and mosaic agriculture.

Figure 1 Diagram of the different LULCC driver datasets that are used as input for JULES-ES simulation.

2.1.2 Integration of multiple datasets (FAO, HYDE, LUH2, and Mapbiomas)

The LUH2 dataset employs an accounting-based method to provide global annual gridded land use states and transitions for the historical period from 850 to 2023 (Chini et al., 2021; Hurtt et al., 2019, 2020). LUH2 considers each grid cell to be naturally forested or non-forested. This potential vegetation area is adjusted based on information about land use activities, such as agricultural expansion and abandonment, wood harvesting, urban development, and shifting cultivation.

LUH2 received the input of the land fractions in agriculture from the HYDE database (Goldewijk et al., 2017). This is a land use dataset that combines historical population estimates and allocation algorithms with time-dependent weighting maps to ensure internal consistency (Goldewijk et al., 2017). The cropland and grazing data in HYDE are obtained from the country-level statistical data provided by FAO as shown in Fig. 1. The cropland in HYDE is from FAO arable land and permanent land, while the grazing in HYDE is from FAO permanent meadows and pastures.

FAO utilized a national statistics approach provided by each country (Ricciardi et al., 2018). The national data are provided to FAO annually via questionnaires completed by each country; however, they may be based on less frequent census periods, which can vary among countries. In Indonesia, FAO has sent out questionnaires since the 1960s. It revised the land use questionnaire in 2016 (FAO, 2022). The data before that year are from the Indonesia Statistic Agency (BPS), which conducted agriculture surveys every 10 years (Badan Pusat Statistik, 2023). The years in between are interpolated. To spatially allocate the FAO data, HYDE uses maps from the ESA Land Cover Consortium (Fig. 1).

Table 1 The LULCC dataset.

This version of LUH2 with cropland areas based on HYDE or FAO (except for Brazil, which started to use Mapbiomas in 2022) was utilized by the Global Carbon Budget until the 2022 version (GCB22) (Chini and Hurtt, 2025). As the FAO dataset is only available until 2020 for GCB22, the years 2021 and 2022 are extrapolated from HYDE using the trend over the previous 5 years (2016–2020). Then, specifically for Indonesia, in the Global Carbon Budget 2023 (GCB23), FAO data for Indonesia were replaced (while other countries except for Brazil continued to be based on HYDE or FAO) with the satellite-based MB1 dataset as input for HYDE and then used by LUH2 for the period 2000–2019, with extrapolation from HYDE until 2022. Recently, Mapbiomas launched the new version, MB2, which is claimed to have higher accuracy. This dataset is currently being proposed for the future GCB. Thus, in this research we differentiate between LUH2-GCB2022, MB1, and MB2. The details of all of the LULCC datasets assessed in this research are given in Table 1.

2.2 ELUC model and product comparison

This study utilizes the JULES-ES model (Sellar et al., 2019). JULES-ES was developed from the Met Office's Surface Exchange Scheme (MOSES). It can be used as a stand-alone land surface model driven by observed forcing data or coupled to an atmospheric global circulation model, such as the Met Office's Unified Model (UM) (Best et al., 2011). As JULES-ES is a process-based model that simulates the fluxes of carbon, water, energy, and momentum between the land surface and the atmosphere (Clark et al., 2011), it can used to analyse the important role of the land surface in the functioning of the Earth system.

JULES-ES has a detailed representation of land surface processes and includes recent developments in surface physical processes (Wiltshire et al., 2020), wood products (Jones et al., 2011), fire representation (Burton et al., 2019), dynamic vegetation (Cox, 2001; Harper et al., 2018), land use and the nitrogen cycle (Wiltshire et al., 2021), plant physiology, and plant functional types. It can simulate the historical evolution of the land carbon cycle under increasing atmospheric CO₂ concentrations and climate change (Harper et al., 2016, 2018). Moreover, it can simulate the effect of land use change when forced with external LULCC datasets. The model description of JULES-ES regarding the impact of land use change on the carbon cycle can be seen in Supplement Sect. S1 (Eqs. S1, S2, and S3). For this study, we use the same setting of JULES-ES as in TRENDY v12 for GCB23 but vary the fraction of agriculture with several datasets as presented in Table 1. The detailed PFT setting of JULES-ES (Hurtt et al., 2019; Chini and Hurtt, 2025) can be checked in Supplement Table S1.

Apart from looking at the effects of various LULCC datasets on ELUC estimates, we also investigate the ELUC uncertainties using different DGVMs. We analyse 18 DGVMs from TRENDY v12 and GCB23, all driven with the LULCC forcing MB1, which has been extrapolated and adjusted with HYDE and LUH2.

In addition to the DGVMs, we use ELUC estimates for Indonesia from three BKMs that contributed to GCB23: H&C (Houghton and Castanho, 2023; Houghton and Nassikas, 2017), OSCAR (Gasser et al., 2020), and BLUE (Hansis et al., 2015). Although the DGVMs and the BKMs BLUE and OSCAR used the LULCC forcing MB1, this results in a different estimation. Furthermore, these two model types have different approaches and parameterizations. BKM methods track the variations in carbon storage within vegetation, soils, and wood products both pre-LULCC and post-LULCC by employing established growth and carbon decay rates over time. Unlike DGVMs, BKMs do not account for the influence of changing environmental factors on vegetation growth rates. Instead of simulating carbon stocks, BKMs rely on direct observational data for carbon densities, which are derived from biome-level values found in the literature and inventory data (Houghton et al., 1983).

We also compare the ELUC estimates from DGVMs and BKMs with the NGHGI of Indonesia, a comprehensive report on GHG emissions and removals compiled by the Indonesian government. The Indonesian NGHGI follows the IPCC (Intergovernmental Panel on Climate Change) definition of land use emissions, which differs from the ELUC definition used by the GCB, as the NGHGI uses terrestrial fluxes occurring on all land that countries define as “managed”. We specifically use data from NGHGI Biennial Update Reports 2 and 3 (BUR2 and BUR3) (Anwar et al., 2021), which cover the time periods relevant for our study, i.e. 2000–2016 and 2000–2019, respectively. Then, we also investigate the FAO GHG dataset for comparison (FAO, 2023; Francesco Tubiello, 2020; Tubiello et al., 2021).

BUR is a periodic national report that follows UNFCCC (United Nations Framework Convention on Climate Change) guidelines for non-Annex-I countries. The GHG inventory quantifies LULCC emissions using Tier-1 and Tier-2 methods according to the 2006 IPCC Reporting Guidelines: emissions are calculated by multiplying the areas of land cover change by specific emission factors for each land cover conversion. The data collections and calculations are done by ministries, bureaus, or agencies. The data are collected by each relevant government organization and then compiled by the Ministry of Environment and Forestry (Anwar et al., 2021). Indonesia released BUR2 in 2018 (Agung Sugardiman et al., 2018) and BUR3 in 2021 (Anwar et al., 2021). BUR3 uses an improved land cover map by reinterpreting land use based on Landsat imagery and revision of relevant sector data. This creates different estimates for LULCC emissions between BUR2 and BUR3 (Anwar et al., 2021).

The FAO GHG inventory summarizes the greenhouse gas emission dissemination. It follows the Tier-1 methods of the IPCC Guidelines for National GHG Inventories (FAO, 2022, 2023). We use the total CO₂ emissions from land use change for comparison with models and other products. Then, CO₂ emissions from “Drained organic soils” and “Fire in organic soils” are added to represent the peat fire emissions and peat drainage. We convert the data from kilotons CO₂ into pentagrams carbon (PgC).

2.3 Analysis

2.3.1 Annual and spatial LULCC

Quantifying the difference in products of LULCC is the first step of the analysis in this research. We compare the temporal changes in LULCC from the datasets LUH2-GCB2022, MB1, and MB2. MB1 classifies agriculture as a combination of oil palm and other agriculture. MB2 has two additional classifications, i.e. rice paddy and pulpwood plantation. The FAO classifies cropland as arable land and land under permanent crops (including oil palms). In the following, we refer to the Mapbiomas category “agriculture” as “cropland” to simplify the comparison across datasets.

The spatial maps of cropland areas are plotted and analysed in order to assess the differences between LULCC datasets. We plot LUH2-GCB2022, whose cropland and grazing land areas are from HYDE based on FAO country statistics. The Mapbiomas dataset is satellite-based; thus, it is already spatially distributed. However, the MB1 and MB2 datasets need to be processed as follows: first reclassify – palm oil, rice paddy, pulpwood plantation, and other agriculture are defined as 100 % crop, while other classifications (forest, non-vegetated areas, and water bodies) are considered 0 % crops. The spatial resolution of Mapbiomas is 30 m, while LUH2 (used as GCB input) has a spatial resolution of 0.5°. We thus regridded the land cover fractions of MB1 and MB2 to 0.5° resolution. This approach has a small impact on the total cropland area, as illustrated in Fig. S1 in Supplement Sect. S2.

2.3.2 ELUC simulations

The LUH2-GCB2022, MB1, and MB2 datasets are used as input to JULES-ES, which is run at 0.5° resolution (Brasika, 2025a, b). This resolution allows us to capture the geographical conditions of Indonesia, which is an archipelago consisting of many islands with considerable topographic variability. To calculate ELUC, we perform three sets of simulations with JULES-ES using cropland areas from the three LULCC datasets. The simulations are performed according to the TRENDY v12 configuration, using a time-varying climate, CO₂, and nitrogen deposition following the TRENDY protocol (Sitch et al., 2024). The first simulation uses LUH2-GCB2022 over the period 1700–2022. The second and third simulations use LUH2-GCB2022 from 1700–1999 and Mapbiomas (MB1 and MB2) that from 2000–2022. JULES-ES is spun up to steady-state conditions by running 50 spin-up cycles of 20 years starting in 1700, and then each set of simulations follows two scenarios, S2 and S3. Scenario S2 uses time-varying CO₂ and climate forcing but constant present-day land use (from the year 2000). Scenario S3 uses time-varying CO₂, climate forcing, and land use varying from 1700 to 2022 (but uses three different datasets (LUH2-GCB2022, MB1, and MB2) from 2000). ELUC is calculated as the difference in net biome productivity at a grid cell level between these two simulations. The ELUC values from these simulations are compared to understand how differences in LULCC from LUH2-GCB2022, MB1, and MB2 impact the emission of carbon.

We also plot annual time series of all the available results from DGVMs, BKMs, and the NGHGI. We use the ensemble mean of TRENDY v12 for the DGVMs, the model BLUE for BKMs, NGHGI data from BUR3, and the FAO GHG inventory. We plot the net carbon emissions to analyse their trends, magnitudes, and patterns. Then, we re-plot the ensemble DGVMs, BKMs (BLUE, H&C, and OSCAR) and NGHGI (BUR2, BUR3, and FAO GHG) by including the peat fire emissions and peat drainage emissions. The NGHGI reports data on peat fire emissions and peat drainage emissions, whereas for DGVMs and BKMs peat fire emissions and peat drainage emissions are added from an external dataset. Peat fire emissions are taken from the Global Fire Emission Database (GFED4s; van der Werf et al., 2017), while peat drainage emissions are a combination of three spatially explicit datasets (FAO peat drainage emissions for 1990–2020, Conchedda and Tubiello, 2020; peat drainage emissions for 1700–2010 from simulations with the DGVM ORCHIDEE-PEAT, Qiu et al., 2021; and peat drainage emissions for 1701–2021 from simulations with the DGVM LPX-Bern v1.5, Lienert and Joos, 2018). We investigate their ELUC trend with linear regressions, magnitude, and variability with the mean and standard deviation, under conditions with or without peat emissions.

3 Results

3.1 Annual land use changes in Indonesia

In this section, we analyse the cropland dynamics over Indonesia for the last 2 decades. There are three datasets, i.e. LUH2-GCB2022, MB1, and MB2. All of the datasets show that cropland increased between 2000 and 2022, but the total change varied between 10 and 20 Mha (Fig. 2a). The highest increase, in LUH2-GCB2022, was mainly in 2016–2022, when it grew by 10 Mha. The difference between the datasets can be seen clearly from the dynamics of the cropland change (Fig. 2b). LUH2-GCB2022 shows higher cropland changes than MB1 or MB2. In certain years, it shows a gain of nearly 3 Mha, such as in 2003 and 2016, while it shows a significant loss of over 1 Mha in 2005, for which we could not find any independent evidence supporting such a large decline in cropland during that year. Furthermore, the order 2010, 2012, and 2014, when the change in LUH2-GCB2022 cropland is exactly zero, seems unlikely as Indonesia has experienced continuous increased development in agriculture during recent decades (Gandharum et al., 2022; Maheng et al., 2021). Other than that, MB1 shows a decrease in cropland area in 2019, but this year has low confidence for MB1, as 2019 is the final year of the dataset and the Mapbiomas methodology needs validation from 1 year before and 1 year after (Mapbiomas Indonesia, 2022). This is confirmed in MB2, which shows no decrease in cropland in 2019. Annual cropland changes in MB2 increased from 2000 to 2011/2012, with a peak of about 1.2 Mha cropland expansion and a slow decrease thereafter.

Figure 2 The comparison of data driver time series of (a) cropland areas (Mha) and (b) annual cropland area change (Mha) as well as the spatial distributions of cropland fractions in (c) LUH2-GCB2022 in 2000, (d) LUH2-GCB2022 in 2018, (e) LUH2-GCB2022 between 2000 and 2018 (2018 fraction minus the 2000 fraction), (f) MB1 in 2000, (g) MB1 in 2018, (h) MB1 between 2000 and 2018 (2018 fraction minus the 2000 fraction), (i) MB2 in 2000, (j) MB2 in 2018, and (k) MB2 between 2000 and 2018 (2018 fraction minus the 2000 fraction).

The differences in cropland change in Indonesia across the datasets can also be observed from their spatial distribution. All three datasets, LUH2-GCB2022, MB1, and MB2, agree that most cropland areas are concentrated on the three main islands of Sumatra, Kalimantan, and Java, with smaller areas found in Bali, Nusa Tenggara, Sulawesi, Papua, and Maluku (Fig. 2c–k). When examining the spatial distribution of the cropland total, MB1 and MB2 show higher cropland fractions in Sumatra, Kalimantan, and Java compared to the LUH2-GCB2022 dataset in 2000 and 2018 (Fig. 2c, d, f, g, i, and j). The LUH2-GCB2022 dataset shows a significantly larger increase in cropland fractions compared to Mapbiomas in Sumatra and Java over the 2 decades (Fig. 2e, h, and k).

However, this does not mean that Mapbiomas is less uncertain than LUH2-GCB2022. For instance, MB1 and MB2 have different spatial distributions to some extent, although they both come from the same source with different versions. For example, in Sumatra in 2000, MB1 shows that a high cropland fraction is distributed across the island (Fig 2f), while MB2 is more concentrated on the southern part of the island (Fig. 2i). This is also seen clearly in Kalimantan (Fig. 2h and k).

3.2 Annual model ELUC estimation

There are two sources of uncertainty in determining the ELUC magnitude and trend in Indonesia: first the driving LULCC datasets and second the fate of the different carbon reservoirs (vegetation, litter, and soils) after LULCC. We first investigate the uncertainty from LULCC by performing simulations with the different cropland datasets (LUH2-GCB2022, MB1, and MB2) in the same DGVM, JULES-ES (Fig. 3a).

Unsurprisingly, JULES-ES driven by LUH2-GCB2022 shows the largest interannual variability in ELUC, consistent with the changes in cropland areas in that forcing (Fig. 2b), showing higher emissions of more than 0.15 PgC yr⁻¹ in 2004 and 2017, and following the largest cropland changes that happened the year before in the LUH2-GCB2022 dataset. The ELUC mean and internal variability from LUH2-GCB2022 is 0.11 ± 0.03 PgC yr⁻¹. These patterns do not appear in JULES-ES when forced by either MB1 or MB2, with the simulated ELUC showing less year-to-year variability and being 0.08 ± 0.01 and 0.06 ± 0.01 PgC yr⁻¹, respectively, for the last 2 decades. However, notable differences also exist between ELUC driven by MB1 and MB2, following the cropland patterns shown in Fig. 2. The MB2 ELUC is lower than the MB1 ELUC, as cropland expansion is generally lower in MB2 (Fig. 2b). The average estimation of the Indonesian ELUC of JULES-ES from different driving LULCC is 0.08 [0.06–0.11] PgC yr⁻¹.

Figure 3 The ELUC model comparison: (a) JULES-ES with various drivers + ensemble in TRENDY v12 and (b) various DGVMs in TRENDY v12 using MB1 as the driver. Spatial and decadal changes of (c) JULES-ES (LUH2-GCB2022) for 2000–2009, (d) JULES-ES (MB1) for 2000–2009, (e) JULES-ES (MB2) for 2000–2009, (f) JULES-ES (LUH2-GCB2022) for 2010–2019, (g) JULES-ES (MB1) for 2010–2019, and (h) JULES-ES (MB2) for 2010–2019.

We now turn to the second source of uncertainty which comes from the representation of ELUC dynamics in models. This is illustrated in Fig. 3b, which shows 18 DGVMs (including JULES-ES) from GCB 2023. All of the DGVMs have the same cropland input MB1, similar to GCB 2023. The result shows that there is little agreement on the magnitude and year-to-year variability of ELUC across the models. The ELUC estimates from various carbon reservoir models in TRENDY v12 are 0.12 ± 0.02 PgC yr⁻¹ (α=0.05). Moreover, there is a large range from −0.04 to 0.35 PgC yr⁻¹, much larger than the range of ELUC estimates from various drivers.

In addition, we can see the impact of different driver data from their spatial distribution. The higher cropland change in LUH2-GCB2022 (Fig. 2e) results in a higher estimation of ELUC, mainly in Kalimantan, Sumatra, and Java (Fig. 3c and f). The JULES-ES simulations based on input from MB1 and MB2 are more similar but still show some differences in the fine spatial details (Fig. 3d, e, g, and h). For example, in Sumatra, JULES-ES (MB2) shows higher carbon emissions in the centre of the island, while the JULES-ES (MB1) simulation has emissions spread more evenly across the island.

3.3 Comparison of ELUC with other products

In addition to analysing the uncertainties caused by land use change drivers and the land cover change dynamics in DGVMs, we compare our findings with other estimates such as those from BKMs and those from the NGHGI. From the Global Carbon Budget 2023, there are three BKMs, which is comparable to the DGVMs analysed above. For the NGHGI, there are BUR2 and BUR3, released by the Ministry of Environment and Forestry of the Republic of Indonesia and covering the 2000–2016 and 2000–2019 periods, respectively. Then the FAO GHG inventory (FAO, 2023; Francesco Tubiello, 2020; Tubiello et al., 2021) is also compared.

The result shows that the ensemble of the TRENDY v12 ELUC mean and internal variability is similar to the bookkeeping model BLUE's ELUC estimation, with 0.12 ± 0.02 and 0.12 ± 0.02 PgC yr⁻¹, respectively, both showing a steady trend (Fig. 4a). The other BKMs H&C and OSCAR result in 0.14 ± 0.05 and 0.15 ± 0.03 PgC yr⁻¹. The BKM mean is 0.14 [0.12 to 0.15] PgC yr⁻¹. The FAO GHG estimates lower carbon emissions of 0.10 ± 0.06 PgC yr⁻¹. The NGHGI's BUR3 estimates the lowest carbon emission mean and internal variability of 0.06 ± 0.06 PgC yr⁻¹, with an increasing trend. Furthermore, BUR3 and the FAO GHG have different annual patterns.

These emissions are much higher when the peat emissions are included. The carbon emissions are doubled: TRENDY v12 = 0.23 ± 0.05, BLUE = 0.22 ± 0.06, and BUR3 = 0.13 ± 0.10 PgC yr⁻¹. This is consistent with other BKMs such as H&C and OSCAR with estimations of 0.25 ± 0.07 and 0.25 ± 0.06 PgC yr⁻¹. The FAO GHG with peat emissions estimates 0.20 ± 0.09 PgC yr⁻¹. Also, another NGHGI, BUR2, estimates 0.18 ± 0.07 PgC yr⁻¹ of the carbon emission total with peat. A distinct difference is also shown by all DGVMs and BKMs having a steady trend, while the NGHGI tends to have an increasing trend and higher variability. We also found that all of the models and products show the same peaks in the years 2002, 2006, 2009, 2015, and 2019 (Fig. 4b). This does not appear in the emissions without peat (Fig. 4a), except for the FAO GHG inventory. This is the exact same year of moderate to high El Niño.

Figure 4 ELUC comparison of different models over the period 2000–2022. (a) Annual ELUC from the DGVMs, BKMs, and NGHGI. (b) Annual ELUC from the BLUE and BUR3 models with and without peat emissions + ensemble TRENDY v12.

4 Discussion

This research analyses land cover products that are used to estimate Indonesian ELUC fluxes. The LULCC driving dataset is one of the main sources of uncertainty in ELUC trends in Indonesia (Bastos et al., 2020; Gasser et al., 2020). Thus, we compare LUH2-GCB22, which is derived from a national-statistics-based dataset of FAO and which is used widely for ELUC, including for the Global Carbon Budget until its 2022 version, with the new satellite-based dataset Mapbiomas Indonesia 1.0, which is utilized as a driver dataset in GCB 2023. In addition, we investigate the latest version of Mapbiomas Indonesia 2.0, which has not yet been used in any GCB assessments. Results show that all of the datasets have an increase in cropland area over the period 2000–2022, but the magnitude and year-to-year change show limited agreement across the datasets. This greatly impacts the simulated annual emissions of ELUC using our carbon cycle model, JULES-ES.

Many models and simulations use changes in agricultural areas based on statistical data reported by countries to FAO as a forcing for tree cover loss. They show that FAO agricultural statistics have high uncertainties (Ajaz, 2016; Desiere et al., 2016). However, the FAO dataset has the capability to be used for net deforestation (Tubiello et al., 2021). On the other hand, the use of satellite-based estimates of land cover changes for land use changes has some limitations and uncertainties (Tubiello et al., 2023). This means that both national-statistics-based and satellite-based datasets have their own advantages and disadvantages, which should be treated with caution in order to improve our understanding of carbon emissions from LULCC.

Cropland transition is essential. For example, in southern Sumatra, where the national-statistics-based (LUH2-GCB2022) and satellite-based (MB1 and MB2) datasets show distinct differences, LUH2-GCB2022 shows a much smaller cropland area compared to Mapbiomas in 2000 (Fig. 2c, g, and i). In this region, forest cover shrank rapidly between 1996 and 2000 (Purnomo et al., 2023). However, this was not directly converted into cropland: it became shrubland in 2000 and was then converted into pulpwood plantation and oil palm from 2000 onwards (Gaveau et al., 2007; Purnomo et al., 2023). This process might have caused different classifications between satellite and national statistics. Furthermore, it can be seen in Java, where a national-statistics-based dataset shows a higher increase in cropland change. Meanwhile, Java has been experiencing massive infrastructure development that limits the potential expansion of agriculture in this period (Gandharum et al., 2022; Maheng et al., 2021).

This Mapbiomas difference with LUH2-GCB2022 aligns with other research that shows Kalimantan and Sumatra having the highest rate of deforestation before 2000 (since 1950) (Gaveau et al., 2022; Santoro et al., 2023). Deforestation in Java and Bali already occurred before 1950 (Gaveau et al., 2022; Santoro et al., 2023). This means that all three of the main islands already had large cropland areas before 2000. There are also large differences in the border regions of Kalimantan and Sumatra. These areas are known for their dense forests and have experienced major forest fires during El Niño events in the last 2 decades (2000–2020) (Brasika, 2023; Fanin and van der Werf, 2017).

The satellite-based dataset has more consistency with cropland expansion in Indonesia, which is related to other LULCC like deforestation and urbanization. This aligns with other research (Gaveau et al., 2022) that shows that forest conversion to agriculture (specifically oil palm) increased during the 2000s, peaked around 2009–2012, and steadily declined thereafter. Also, MB2 does not show peculiar years when cropland area changes would be exactly zero. Overall, we expect MB2 to be the most accurate representation of the croplands we assess in this study.

In our current model, we use grazing lands from the LUH2-GCB2022 dataset as input for all simulations in JULES-ES. This might cause some double accounting of grazing land in our satellite-based simulations, where grazing land might be counted as both cropland and grazing land. However, this effect should be very small as grazing land covers a small area in Indonesia (around 3 Mha) and has been considerably steady for the last 2 decades (FAO, 2022; Mapbiomas Indonesia, 2022). Also, in the JULES-ES simulation, there is no large difference in the carbon emissions caused by forest conversion to cropland or grazing land (Burton et al., 2019).

Furthermore, another potential misrepresentation in our model is ELUC from palm oil. As palm oil makes significant contributions to the forest–agricultural dynamic in Indonesia (Cisneros et al., 2021; FAO, 2023; Gaveau et al., 2021; Tsujino et al., 2016), it should be represented in the model as a specific plant functional type (PFT). Unfortunately, we currently follow the Mapbiomas Indonesia definition of agriculture, which categorizes palm oil (together with rice paddy, pulpwood plantation, and other agriculture) as cropland, because there is no PFT for palm oil in our model. The lack of a palm oil vegetation type is likely to cause an inaccuracy of ELUC. As Mapbiomas Indonesia has palm oil in Indonesian land cover (Mapbiomas Indonesia, 2022), future model simulations should consider this separation.

Another interesting aspect is how these products might reveal the connection between carbon emissions and climate dynamics. All DGVMs, BKMs, and the NGHGI with peat emissions show a strong relation to climate events, which shows that the peak in carbon emissions is in exactly the same year of moderate to strong El Niño. This appears in 2002, 2006, 2009, 2015, and 2019. It is mainly contributed by the carbon emissions from peat fire and peat drainage. During the El Niño year, Indonesia experiences drier and hotter climate conditions, especially in Sumatra and Kalimantan (Brasika, 2021; Nurdiati et al., 2022). This condition is favourable for fire regimes in the area that burn not only the above-ground biomass, but also the peat soil below the ground (Brasika et al., 2021; Fanin and van der Werf, 2017). As peat contains a lot of carbon and peat fire is hardly detected and managed (Indradjad et al., 2024), this results in the massive peat fire release of carbon during the hot and dry El Niño year (Stockwell et al., 2016). However, ELUC without peat emissions has almost zero connection to these El Niño dynamics. This is confirmed by all of the products of the DGVMs, BKMs, and NGHGI. The carbon emissions from peat play a key role in Indonesia, with fluxes 1–3 times higher than ELUC and stronger in El Niño years. This should not be neglected in future simulations.

5 Conclusions

Both the national-statistics-based and satellite-based datasets have their own strengths and weaknesses. The national-statistics-based dataset has lower confidence as it has no spatial distribution, while the satellite-based dataset has limited temporal availability. Addressing these differences might improve our understanding of uncertainties in carbon emission estimates caused by the drivers. However, improving the driver dataset is not the only issue: the model itself has potential uncertainties as each of them is created with different approaches. This can be seen in how all of the DGVMs in TRENDY v12 have the same input drivers but result in different patterns and magnitudes of carbon emission simulations. Our best estimation of Indonesian ELUC for the last 2 decades is 0.12 ± 0.02 PgC yr⁻¹. As this is the mean of all of the DGVMs, it is also consistent with BKMs. This emission is doubled to 0.23 ± 0.05 PgC yr⁻¹ if we include the peat fire and peat drainage emissions.

We find some important issues that should be addressed carefully in order to reduce these uncertainties in the future. First, the combination of the national-statistics-based and satellite-based datasets improves our understanding of the uncertainties of carbon emissions. Second, the inclusion of peat emissions and the climate effect (El Niño) in the models makes major contributions to the carbon emissions in Indonesia. Third, palm oil representations, which are dominant in Indonesia's agriculture with a distinct carbon cycle, have been simplified here.

Apart from the many uncertainties created by the drivers and models and some representations such as peat emissions, palm oil, and climate dynamics, we found several robust agreements between these models and products. First, all of the products agree that Indonesian carbon emissions from LULCC have had no decreasing trend for the last 2 decades. Second, the peat fire and peat drainage have doubled the carbon emissions. Third, it is highly confident that carbon emissions of peat have a strong relation to El Niño but not carbon emissions from land use change. Fourth, Sumatra, Kalimantan, and Java are the main islands which contribute to Indonesian carbon emissions. Last, it seems that prevention of deforestation and fire had a low effect in Indonesia in reducing carbon emissions when addressing climate change.