Uncertainty in Land Carbon Fluxes Simulated by CMIP6 Models from Treatments of Crop Distributions and Photosynthetic Pathways

A reliable representation of the diversity of vegetation in terrestrial ecosystems is needed for the accurate simulation of present and future biogeochemical cycling and global climate, particularly as climate change affects different vegetation types differently. We compare the distributions of crops and of C₃ vs. C₄ photosynthetic pathways in both natural vegetation and crops across Earth System Models in the 6th Coupled Model Intercomparison Project (CMIP6). We find a large range in C₃ and C₄ crops, natural and total vegetation for area and gross primary production (GPP) across the models. Even though 10 of the 11 models used Land Use Harmonisation (LUH2) crop areas as input data, modelled total crop area ranges from − 28 % to + 10 % of a satellite-based estimate. The C₃ and C₄ crop areas were − 56 % to + 15 % and − 100 % to + 38 % of LUH2 for 2014, respectively. The C₄ fraction of total vegetation area in the models is 9 %–25 %, compared to 20 % ± 3 % in observation-based estimates for the year 2014. Total global GPP varies by a factor of two across the models, and the C₄ fraction of GPP ranges from 12 % to 27 %. Simulated trends in the fraction of GPP by C₃ vs. C₄ vegetation type (− 20 % to + 29 %) would have changed global stable carbon isotopic discrimination by − 0.35 ‰ to + 0.11 ‰ over 1975–2005, not including changes in discrimination over time within C₃ and C₄ vegetation, indicating that modeled changes in the fraction of GPP by C₃ and C₄ vegetation do not account for the + 0.7 ‰ increase indicated by atmospheric data. Disparity in vegetation with these photosynthetic pathways in models contributes to uncertainty in land carbon flux simulations, and further constraints and improvements in models are needed.

1 Introduction

The terrestrial biosphere captures ∼ 30 % of anthropogenically emitted CO₂ annually (Friedlingstein et al., 2025) reducing CO₂ accumulation in the atmosphere and the accompanying global warming. However, this CO₂ uptake may be sensitive to future climate change (Arora et al., 2020). Understanding the mechanisms contributing to the CO₂ uptake, including the role of different types of vegetation and land use, is essential to understanding potential carbon-climate feedback and future changes to the terrestrial carbon cycle.

The implementation of land use and land cover change (LUCC) and its impact on vegetation cover and dynamics are important components of terrestrial biosphere model development (Hu et al., 2021; Lam et al., 2021; Hurtt et al., 2020; Wang et al., 2022). LUCC is mainly driven by agriculture. The fraction of global land area used for agriculture increased from 14 % in 1850 to about 37 % in 2015 (Hurtt et al., 2020), and LUCC emitted about 118 PgC between 1850 and 2020 (Houghton and Castanho, 2023). In addition to CO₂ emissions, the biophysical effects of land conversion also drive global temperature rise due to changes in surface albedo (Arora and Boer, 2010; Houghton et al., 2012). The alteration of the land surface will continue to be significant in the future as, for example, it has been projected that 14 % of vegetation and 5 % of soil carbon stocks will be lost to cropland expansion globally over 2010–2050 in a “middle-of-the-road” scenario (Molotoks et al., 2018).

To represent plant diversity and function, most model developers use plant functional types (PFTs) that group vegetation by similar features such as growth form, ecological requirements, and photosynthetic pathways. This helps to account for the variation in adaptive mechanisms and ecological distribution of different plants (Haxeltine and Prentice, 1996; Hurtt et al., 2020; Wullschleger et al., 2014). One important characteristic of plants is their use of either the C₃ or C₄ photosynthetic pathways. In C₃ plants, the first product of the photosynthetic pathway is a three-carbon molecule called 3-phosphoglycerate (3-PGA), whereas in C₄ plants, it is a four-carbon molecule called oxaloacetate. C₃ and C₄ types also differ in their response to changes in soil moisture content, temperature, CO₂, and light (Luo et al., 2024). For more information about the anatomical/structural and physiological differences between vegetations with these photosynthetic pathways, which drive their distinct climate and CO₂ concentration responses, see Luo et al. (2024), Cortés et al. (2021), Smith and Boers (2023), Polley et al. (1994), Farquhar et al. (1989), Still et al. (2003), Ehleringer and Björkman (1977).

Under increasing temperatures, especially those above the thermal optimum threshold, productivity in C₃ vegetation is limited by increased photorespiration (Hermida-Carrera et al., 2016), which reduces photosynthetic efficiency in C₃ plants. In contrast, C₄ species overcome photorespiration through their CO₂ concentrating mechanism that increases the amount of CO₂ at the site of carboxylation (da Silva et al., 2020), leading to an abundance of C₄ species (natural grasses and crops such as maize) in hot areas in the tropics and sub-tropics (Luo et al., 2024). While rising temperature generally favours C₄ plants, rising atmospheric CO₂ concentration confers a physiological advantage upon C₃ species (Polley et al., 1994) due to a reduction in photorespiration and an increase in water use efficiency. This CO₂ fertilisation effect is responsible for the increasing presence of C₃ woody species in previously C₄ dominated grasslands (Luo et al., 2024). Shifts in C₃ and C₄ species composition and carbon fluxes across different regions are projected to continue in future due to changing temperature, water availability and increasing CO₂ concentration in the atmosphere (Cortés et al., 2021; Smith and Boers, 2023; Winslow et al., 2003).

A change in the relative contributions of C₃ and C₄ vegetation to global productivity may contribute to a global trend in stable photosynthetic carbon isotope discrimination (Δ) because C₃ plants discriminate against carbon-13 more strongly than C₄ plants (Farquhar et al., 1989). Since atmospheric studies have indicated that Δ increased by 0.7 ‰ over 1975–2005 globally (Keeling et al., 2017), and by 0.4 ‰ over 2000–2011 in the Northern Hemisphere (Peters et al., 2018), understanding the effect of changes in vegetation type on discrimination would help to quantify the environmental and physiological factors influencing plant function and resulting discrimination, including soil moisture content, vapour pressure deficit and stomatal conductance (Cornwell et al., 2018; Griffis et al., 2010; Francey and Farquhar, 1982).

Currently, the relative contribution and its change over time of C₃ vs. C₄ vegetation to global terrestrial biosphere productivity, and their ecological roles in climate change mitigation, are not well known. Using remote sensing products, physiological modelling, and crop data, Still et al. (2003) found that C₃ and C₄ area abundances are 87.4 and 18.8 million km² (17.7 % as C₄), and C₃ and C₄ gross primary production were 114.7 and 35.3 PgC yr⁻¹ (23 % as C₄), respectively, on average for the 1980s and 1990s. More recently, Luo et al. (2024) used global observations of plant photosynthetic pathways, satellite remote sensing, and photosynthetic optimality theory to estimate a similar level of C₄ area coverage (17.5 %), but a lower fraction of gross primary productivity (19.4 %), compared to Still et al. (2003). Luo et al. (2024) showed that C₄ vegetation coverage decreased from 17.7 % to 17.1 % over 2001 to 2019 as natural C₄ grass cover declined in favour of C₃ vegetation, especially C₃ trees in tropical grasslands and savannas. In comparison, across the TRENDY ensemble of dynamic global vegetation models, there were large ranges of 7 %–23 % of vegetation area and 2 %–40 % of productivity from C₄ vegetation (Luo et al., 2024), showing a need for better constraints and understanding of C₃ $/$ C₄ vegetation competition and change over time.

Here we investigate the contribution of vegetation diversity representation to uncertainty in carbon flux simulation in Earth System Models. We evaluate the representation of C₃ and C₄ vegetation over 1850–2014 in 11 Earth System Models used in the 6th Coupled Model Intercomparison project (CMIP6). We assess the distribution, productivity and carbon content of C₃ and C₄ vegetation and compare with observation/satellite-based estimates where possible. We also explore the potential trend in global carbon isotope discrimination Δ¹³C due to C₃ and C₄ vegetation changes simulated in the models.

2 Materials and Methods

2.1 CMIP6 Models Outputs and the LUH v2 Data

To analyze the CMIP6 models, we obtained the output for the following variables from the CMIP6-ESGF repository (https://esgf-node.ipsl.upmc.fr/search/cmip6-ipsl/, last access: 28 July 2026): fractions of C₃ and C₄ vegetation coverage (c3PftFrac and c4PftFrac), C₃ and C₄ crop fraction (cropFracC3 and cropFracC4), total crop fraction (cropFrac), gross primary production (gpp), grid cell area (areacella), and percentage of each grid cell covered by land (sftlf). The necessary output was available from 11 models (Table 1). Brief descriptions of the implementation of vegetation abundance in each of the 11 models are provided in the Supplement (section “Materials and Methods S1”).

Table 1 Summary of each model's definition of crop cover and our calculation of C₃ and C₄ crop cover. More detail is given in SM Materials and Methods.

To specify crop cover, ten of the eleven CMIP6 models use the Land Use Harmonisation version 2 (LUH2) dataset (Hurtt et al., 2020, Table 1), so we also analyse the LUH2 data directly here. LUH2 estimates the fractional area and transition of use of 12 categories at an annual and 0.25° by 0.25° spatiotemporal resolution, starting from the year 850 (Hurtt et al., 2020). Five of these categories are C₃ annual crops, C₃ perennial crops, C₃ nitrogen-fixing crops, C₄ annual crops and C₄ perennial crops. We group these categories into C₃ crops and C₄ crops, and their sum as total crops.

The UKESM1 model uses the LUH2 crop cover data with some modifications, including simulation of competition between C₃ and C₄ vegetation using the dynamic vegetation model TRIFFID (Top-down Representation of Interactive Foliage and Floral Including Dynamics, Sellar et al., 2019, Clark et al., 2011). Other models use the LUH2 data more directly, but differences can result from mapping the LUH2 data onto the model's PFTs, land cover and grid maps.

In CNRM-CM6.1, the crop cover is based on the ECOCLIMAP-II database (Voldoire et al., 2019). CNRM-CM6.1 adopts a vegetation cover map that is fixed to its present-day distribution, so that temporal changes in crop area are not included. The reasoning for using a fixed distribution is that the carbon cycle is not fully resolved in the model and there is uncertainty in the impacts of land use and land cover change on carbon flux (Faroux et al., 2013; Voldoire et al., 2019; Gasser et al., 2020).

2.2 Global Land Vegetation Cover and Productivity

We analysed changes in the vegetation cover (C₃, C₄ cropland and natural vegetation) for the 11 CMIP6 models by calculating the gridded annual area fraction for C₃ and C₄ vegetation for each year 1850–2014, which we also separated into crops, natural and total vegetation. We also calculated gridded annual gross primary production (GPP) for C₃ and C₄ vegetation in crops, natural and total vegetation in each model. We combined gridded estimates into global totals.

We calculated annual maps of the area fraction of C₃ crops (cropFracC3) and C₄ crops (cropFracC4). Most models included the cropFracC3 and cropFracC4 variables that could be used directly (Table 1). Some models had to be treated differently due to the output provided. ACCESS-ESM1.5, CMCC-CM2-SR5, and CMCC-ESM2 did not provide cropFracC4, and their cropFracC3 was equal to the total crop fraction (cropFrac), so we specified the fraction of C₄ crops to be zero everywhere. The CESM2 models did not provide cropFracC3 and cropFracC4 variables, but they did provide c3PftFrac or c4PftFrac (the fractions for total vegetation), so cropFracC3 and cropFracC4 were calculated by multiplying cropFrac by c3PftFrac or c4PftFrac. In CNRM-CM6-1 and CNRM-ESM2.1, the cropFracC3 and cropFracC4 variables were identical to the c3PftFrac and c4PftFrac variables, which appears to be an error in the output. So, in CNRM-CM6-1 and CNRM-ESM2.1, cropFracC3 and cropFracC4 were calculated by multiplying cropFrac by c3PftFrac or c4PftFrac.We also compared the total cropland area in all the models and the LUH2 dataset to a satellite-based estimate from Potapov et al. (2022).

To calculate the fractions of natural C₃ and C₄ vegetation, we subtracted cropFracC3 and cropFracC4 from c3PftFrac or c4PftFrac (the total C₃ and C₄ vegetation fractions). Therefore, for each grid cell, we had the fractional area of C₃ and C₄ crops and of C₃ and C₄ natural vegetation.

The GPP from C₃ or C₄ vegetation was not provided, so we estimated these by multiplying the total GPP by the area fractions. In Luo et al. (2024), it was found that the per-unit-area photosynthetic rate of C₄ grass was generally higher than that of C₃ vegetation, indicating that these calculations may overestimate the proportion assigned to C₃ vegetation. Therefore, we conducted a reliability analysis to determine the similarity between the values and distribution of C₃ and C₄ crops obtained from cropFrac × c3PftFrac and cropFrac × c4PftFrac in the models (CESM2, CESM2-WACCM, CNRM-ESM2-1 and CNRM-CM6.1) and LUH2 values and distribution using kernel density estimate (KDE) (Silverman, 1986; Chen, 2017) and constructed spatial probability distribution plots using all land pixels for C₃, C₄ and total crops for the year 2014 (Chiang et al., 2021). We also visually compared the spatial distribution of these values to LUH2 and to UKESM1 and MPI-ESM-1-2-HAM, two models that provided their C₃ and C₄ crops.

2.3 Global Stable Carbon Isotopic Discrimination

To estimate how the simulated changes in C₃ and C₄ fractions of GPP influenced the global stable carbon isotopic discrimination, we estimated discrimination assuming fixed values for stable carbon isotope discrimination by C₃ and C₄ vegetation. This neglects any environmental or physiological effects on discrimination that may be simulated by the models to isolate the potential effect from changes in C₃ and C₄ fractions of GPP alone. The total annual global stable carbon isotopic discrimination (Δ_tot) was calculated by the following weighted average with GPP:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{\Delta }}_{\mathrm{tot}}& =({\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{3}}}{\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{3}}}+{\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{4}}}{\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{4}}})\\ & /({\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{3}}}+{\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{4}}})\end{array}\end{array}$$

Here, GPP${}_{{\mathrm{C}}_{\mathrm{3}}}$ and GPP${}_{{\mathrm{C}}_{\mathrm{4}}}$ are the integrated GPP for all C₃ and C₄ vegetation each year, and ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{3}}}$ and ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{4}}}$ are the fixed stable carbon isotope discrimination for C₃ and C₄, respectively. We also calculated the global natural stable carbon isotopic discrimination (Δ_nat) using the natural C₃ and C₄ GPP for all the models using the same weighted average approach. To calculate the effect of crops on the global annual stable carbon isotope discrimination (Δ_crop), we subtracted Δ_nat from Δ_tot.

${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{3}}}$ and ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{4}}}$ were specified as the means of C₃ and C₄ plant leaf stable carbon isotope discrimination observations from the database of 3987 species published by Cornwell et al. (2018), which are 20.7 ‰ for ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{3}}}$ and 6.3 ‰ for ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{4}}}$ (Fig. 1). The two peaks in the histogram are assumed to correspond to the highest count for the species with either the C₃ and C₄ photosynthetic pathways. A breakpoint of 12.5 ‰ was chosen for C₃ and C₄ based on the midpoint between the two peaks, such that discrimination values lower than 12.5 ‰ in Fig. 1 are taken to correspond to vegetation with the C₄ photosynthetic pathway while values higher than 12.5 ‰ correspond to stable carbon isotope discrimination in vegetation with the C₃ photosynthetic pathway. The means of these two distributions were calculated and used to specify the fixed C₃ and C₄ discrimination: 20.7 ‰ for ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{3}}}$ and 6.3 ‰ for ${\mathrm{\Delta }}_{{\mathrm{C}}_{\mathrm{4}}}$.

Figure 1 Stable carbon isotope discrimination in leaves from the database published by Cornwell et al. (2018).

3 Result

3.1 Crops Spatial Distribution and its Change between 1970 and 2014

Based on the LUH2 dataset averaged over 1970–2014, the highest C₃ crop abundance is in agricultural regions in central North America, southwest Europe, southeast and south Asia, Sub-Saharan Africa, southeast South America and southern Australia (Figs. 2 and 3) (Hurtt et al., 2020). The fractional crop area is between 60 % and 100 % in these regions. For C₄ crops, the highest abundance is in the Great Plains in North America, Sub-Saharan Africa, Southwest Asia and Southeast Asia (Fig. 2).

The models that do not use LUH2 data directly, UKESM1 (which used LUH2 indirectly) and CNRM-CM6-1, show differences with LUH2 (Figs. 2 and 3). UKESM1 underestimates C₃ crops in Sub-Saharan and central Africa and in Southeast Asia, particularly in India, which is likely due to decreased precipitation in UKESM1 in these areas since the vegetation cover is coupled to its simulated climate (Sellar et al., 2019) and UKESM1 underestimates precipitation in India. UKESM1 estimates more C₄ crop coverage than LUH2, particularly in Asia (Figs. 2 and 3 and Table S1 in the Supplement). Crop coverage in UKESM1 is about 10 % lower than LUH2 between the equator and 30° N (Fig. 3). CNRM-CM6-1 has less C₄ crop cover in Europe and Africa, but more in North and South America, compared to LUH2 (Figs. 2 and 3). For models using LUH2 directly, there was still disagreement between the models and the LUH2 crop coverage, which could be due to a number of reasons, such as preprocessing of the LUH2 vegetation classes into the plant functional types (PFTs) and remapping of the LUH2 data into the spatial resolution needed for each model (Table S1). For the models where we calculated C₃ and C₄ crop cover by multiplying cropFract by C3PftFrac and C4PftFrac, CESM2, CESM2-WACCM, CNRM-ESM2-1 and CNRM-CM6.1, we found that the kernel density estimate (KDE) and visual comparison of the values of C₃ and C₄ crops distribution show they are within the range of values for UKESM1, MPI-ESM-1-2-HAM and LUH2 in terms of magnitude, range and spatial distribution as shown by the spatial distribution, and the overlap area in the percentage cover vs. density plots. For example, the overlap area between LUH2 and CESM2 in C₃ and C₄ crop cover is 0.68 and 0.70, respectively, which is close enough to that between LUH2 and MPI-ESM-1-2-HAM, which is 0.74 and 0.50, respectively (Figs. 3, S7–S13).

The models incorporating LUH2 have similar patterns in C₃ and C₄ crop; however, the total area of crops in these models differed from LUH2 by − 19 % to + 3 % and from satellite-based estimate by − 28 % to + 10 % (Fig. 4, Table S1) from Potapov et al. (2022) study that analysed the change in global cropland area in the first two decades of the twenty-first century using satellite and sample data that estimate crop cover at ∼ 10.9 and ∼ 11.42 million km² respectively between 2003 and 2006. CanESM5 was the most consistent with LUH2. Inconsistencies among the models that use the LUH2 dataset may be due to how the data were pre-processed before incorporation into the model, due to grid spacing or due to differences in the model PFTs compared to LUH2 categories, particularly as some models did not include C₄ crops. Also, because the variables cropFracC3 and cropFracC4 were not available for all models, our crop-fraction calculations may have contributed to the discrepancies. The magnitude of area abundance for C₃ and C₄ crop and natural vegetation for the years 1970 and 2014 are provided in the Supplement (Table S2).

Between 1970 and 2014, there were decreases in crop coverage in North America, Europe, Southern Africa, Chile, Japan and New Zealand, but increases in South and Southeast Asia, Africa and South America (Figs. 3b and S1). The largest increases of 10 % or more were concentrated around the eastern region of South America, especially in Brazil, the southern fringes of the Sahara Desert (sub-Saharan Africa) and southeast Africa. The largest decreases of more than 20 % were across the corn belt of the United States and in southwest and central Europe. Changes in crop area in UKESM1 were consistent with LUH2 except in the low latitudes of the Northern Hemisphere, where increases in C₃ crops in Africa and Asia were underestimated, compared to LUH2 (Figs. 3b and S1).

Figure 2 Mean percentage of Earth surface covered by crops between 1970 and 2014 for C₃ and C₄ in LUH2, CNRM-CM6.1 and UKESM1-0-LL. Other models not shown use the LUH2 crop variables directly for their crop coverage and change.

Figure 3 Mean and change in crop coverage averaged over latitude between 1970 and 2014. (a) Mean crop coverage (b) change in crop coverage.

3.2 Global Temporal Trends in Vegetation Coverage

There is a strong positive trend in the area of croplands from 1850 until 2014 (Fig. 4a). Total crop area rose by 200 % in LUH2 and models using LUH2 data. Crop area in UKESM1 was lower than LUH2 but increased in a greater proportion, from 3 to 10 million m² (217 %). In CNRM-CM6, the total crop area was fixed at ∼ 14 million km², similar to LUH2 in the 1990s.

While representing a smaller fraction of total crop area (10 %–45 %) (Fig. 4a and g), the area of C₄ crops increased in greater proportion (> 300 %) than C₃ crops (160 %) in the LUH2 data (Fig. 4). Historically, in the CESM2 models, C₃ crop area is lower and C₄ crop area is higher than in LUH2, despite the model using LUH2 data, although this may be affected by our calculation from the variables provided (Table 1). In UKESM1, the C₃ crop area was much smaller than LUH2, while its C₄ crop area was larger (Fig. 4). UKESM1 C₃ crop area peaked between the 1980s and 90s before declining slightly in the 2000s, in contrast to LUH2, where it continued to increase after 2010. UKESM1 C₄ crop area rose through 2014, when it was about 42 % higher than LUH2. CNRM-CM6-1 had a fixed C₄ crop fraction of 24 %, slightly higher than the C₄ crop fraction in LUH2 in the 1990s. ACCESS-ESM1-5, CMCC-ESM2 and CMCC-CM2-SR5 do not have C₄ crops, so their total crop area is equivalent to the C₃ crop area.

The area of natural vegetation decreases in all models from 1850 until 1940–1970, when four models (CNRM-ESM2-1, UKESM1, MPI-ESM-1-2-HAM and MPI-ESM1-2-LR) start increasing while the others continue decreasing. Natural vegetation is similarly dominated by C₃, with 9 %–23 % of natural vegetation as C₄ in 2014, and most of the decline in natural vegetation area is in C₃ vegetation. The trend in natural C₄ vegetation area is inconsistent across the models. In CanESM5 and ACCESS-ESM1-5, the natural C₄ vegetation area decreased over 1850–2014, while for MPI-ESM-1-2-HAM and MPI-ESM1-2-LR the natural C₄ vegetation area increased especially from 1990. For the CESM2, CMCC-ESM2 and CMCC-CM2-SR5 , the natural C₄ vegetation area was constant until the early 2000s before falling slightly. The UKESM1 estimate of the area of natural C₄ vegetation for 2000–2014 compares well with Luo et al. (2024)'s, while other models simulate up to 15 % larger and up to 63 % smaller areas (Fig. 4h).

Figure 4 Global temporal trends in vegetation area from 1850 to 2014 in CMIP6 models. Trend of the area covered by (a) total crop, (b) total natural vegetation (c) total global vegetation, (d) C₃ crops, (e) C₃ natural vegetation, (f) total C₃ vegetation, (g) C₄ crop, (h) C₄ natural vegetation and (i) total C₄ vegetation. All models except CNRM-CM6-1 used LUH2 data to inform their crop coverage. ACCESS-ESM1-5, CMCC-ESM2 and CMCC-CM2-SR5 do not have C₄ crops.

For the total vegetation area, four models have distinct positive trends from the late 1960s through to 2014 (CNRM-ESM2.1, UKESM1, MPI-ESM-1-2-HAM and MPI-ESM1-2-LR) (Fig. 4c), driven mostly by an increase in C₄ natural vegetation (Fig. 4h). In the other models, the increase in crop area is balanced by the decrease in area covered by natural vegetation. The estimated total vegetated area from the European Space Agency Climate Change Initiative (ESA-CCI) for 2000–2014 is matched by the CMCC-ESM2 and CMCC-CM2-SR5, while the CESM2 models simulate larger vegetated areas and all other models simulate smaller vegetated areas. The ESA CCI data shows a small increase that may be caused by the replacement of bare ground by natural grasses (Fig. S2).

3.3 Global Trends in Gross Primary Production

While all the models simulate an increase in GPP, the magnitude and contribution by vegetation types (in this context: crop/natural and C₃ $/$ C₄) differ (Fig. 5). The models disagree on the magnitude of total global GPP, ranging from 80 PgC yr⁻¹ in CNRM- ESM2.1 to 150 PgC yr⁻¹ in MPI-ESM-1-2-HAM and MPI-ESM1-2-LR (− 20 % to + 29 %) in 2014 (Gier et al., 2024; Arora et al., 2020) (Fig 5c). The total GPP for all the models increased steadily until the 1960s and then grew sharply for the rest of the historical period (Fig. 5c).

The increase in GPP in crops is linked to the cropland expansion, while the GPP change in natural vegetation is decoupled from the change in area (Figs. S5 and S6). In Fig. 4a, the total crop area increased by 200 % compared to Fig. 4b, in which natural vegetation decreased by 7 % on average. Comparing this change in area to the change in GPP in Fig. 5a and b, total crop GPP increased by 209 % on average. However, for CNRM-CM6.1 with a temporally fixed vegetation cover, the increase in total crop GPP is only 18 %. For natural vegetation, even though there is a mean 7 % decrease in total natural vegetation cover, there is an ∼ 18 % increase in GPP on average. In CNRM-CM6-1 with a fixed vegetation cover, the trend of C₃, C₄ and total crop GPP is not as strong compared to other models, especially before 1970 (Fig. 5a, d and g). They increased by 32 %, 27 %, and 18 % for C₃, C₄ and total crop GPP, respectively, compared to the ensemble mean increase of 195 %, 251 % and 209 %, respectively (Fig. 5a, d and g).

For GPP in natural vegetation, there are large ranges of 63–107 PgC yr⁻¹ simulated for C₃ vegetation and 9–23 PgC yr⁻¹ simulated for C₄ vegetation across the models before 1970 (Fig. 5e and h). Despite the decrease in natural C₃ vegetation area before 1970, the C₃ GPP trend is generally either weakly negative or unchanging until 1970 (Fig. 5e). Then all models increase steadily over 1970–2014. In natural C₄ vegetation, there is generally an increase in GPP since 1850, but the increase is not at the same rate in all the models. It is weaker in ACCESS-ESM1-5, CanESM5 and UKESM1 (Fig. 5h). Natural C₄ GPP decreased and remained relatively constant in ACCESS-ESM1-5 and CESM2 models, respectively, over 2000–2014.

The proportion of total GPP by C₄ vegetation in the models is 12 %–27 % (Fig. 5i), compared to 23 % in Still et al. (2003) and 20 ± 3 % in Luo et al. (2024). Compared to other model simulations, the range was 2 %–40 % in TRENDY models (Luo et al., 2024), and 18 %–27 % in previous modelling studies (Farquhar and Lloyd, 1993; Fung et al., 1997).

Figure 5 Temporal trend in gross primary production (GPP) between 1850 and 2014. Trend of GPP of (a) total crop, (b) total natural vegetation (c) total global vegetation, (d) C₃ crops, (e) C₃ natural vegetation, (f) total C₃ vegetation, (g) C₄ crop, (h) C₄ natural vegetation and (i) total C₄ vegetation.

3.4 Global Trends in Stable Carbon Isotope Discrimination

The models disagree on the magnitude and trend of Δ_tot, Δ_nat and Δ_crop (Fig. 6), for our calculation of discrimination based only on the fraction of GPP from C₃ and C₄ vegetation, neglecting any physiological or environmental effects Eq. (1). As a result of the 12 %–27 % range in C₄ fraction of total GPP, the global total stable carbon discrimination ranges from 16.9 ‰ to 19.5 ‰ (Table S2) in 2014 across the models, compared to 16.5 ‰ in Still et al. (2003). Models simulate no change (CESM2 models), slightly increasing Δ_tot (ACCESS-ESM1.5 and CNRM-CM6.1) or decreasing Δ_tot (all other models). The strongest decreases in Δ_tot were in MPI-ESM-1-2-HAM and MPI-ESM1-2-LR, which had the strongest increases in C₄ vegetation area (Fig. 6). Compared to the strong positive trend in the global stable carbon isotope discrimination trend derived by Keeling et al. (2017) from atmospheric δ¹³C data, the effect of vegetation change in the models was much smaller or had the opposite sign (Fig. 6).

The models also disagree on the sign and magnitude of the trends in change in Δ_nat (− 0.83 ‰ to 1.7 ‰) and Δ_crop (− 0.5 ‰ to 0.2 ‰ , Table S2). Among models with decreasing Δ_tot, the trend is primarily driven by growth in the productivity of C₄ natural vegetation, affecting Δ_nat in MPI-ESM-1-2-HAM and MPI-ESM1-2-LR, but by growth in the productivity of C₄ Crops, affecting Δ_crop in the UKESM1 and CanESM5 models (Fig. 6).

Figure 6 Trend of global stable carbon isotope discrimination between 1850 and 2014. 10 year average stable carbon isotope discrimination of (a) total vegetation (Δ_tot), (b) difference in Δ_tot from 1850, (c) component of difference due to natural vegetation (Δ_nat) and (d) due to crops (Δ_crop).

4 Discussion

By examining the components of global vegetation area and GPP by vegetation types, we find that CMIP6 models span a large range for nearly all variables. Some models showed stronger gains in GPP in natural C₃ vegetation, while for others the strongest gains were in either crop or natural C₄ vegetation. Therefore, better quantification of vegetation types and more consistency in models would improve future carbon modelling. For example, total vegetated areas in some models that are much lower than the ESA-CCI estimate (Harper et al., 2023) and total C₄ vegetated areas in some models that are much lower than the Luo et al. (2024) estimate can probably be ruled out. Additionally, Zhao et al. (2025) show that the uncertainty in the terrestrial biosphere GPP simulation in models depends strongly on the area of each PFTs and this dependence enabled them to reduce the spread in land carbon estimates by ∼ 75 %. They also found that changes in PFTs distribution is responsible for 56 % ± 21 % climate-induced variability in dynamic global vegetation models (DGVMs) GPP. This spread across models in many cases reflects the lack of observational constraints and may also be linked to the discrepancies in the natural vegetation composition in the context of our study.

Even though nearly all the models use LUH2 as input data, their crop areas and, in particular, their attribution to C₃ and C₄ vegetation are not consistent. While most models underestimated total crop area compared to LUH2 (Fig. 4a), another satellite-based cropland map (Potapov et al., 2022) estimates a lower crop area that overlaps some of the CMIP6 models' crop areas. Therefore, the spread in crop area across the models may reflect uncertainty in their input crop-area data. A limitation of some current models is that they do not include C₄ crops; thus, carbon fluxes associated with C₄ crops cannot be simulated in those models. In the future, the total crop area is likely to increase (Molotoks et al., 2018; O'Neill et al., 2016), so agriculture will exert an even stronger influence on carbon fluxes.

The range in global GPP was even larger than the range in vegetated area across the models, both for total (GPP: − 23 % to + 37 % vs. area: − 11 % to + 8 %, compared to the model mean) and for C₃ (GPP: − 24 % to + 32 % vs. area: − 19% to + 8 %) vegetation types. In C₄ natural vegetation, there was also disagreement over the sign of the GPP trend, with some models showing increases in GPP in grid cells dominated by C₄ (> 75 %; Fig. S3; a–d), and others showing decreases. These inconsistencies indicate that models' parametrisations, simulated climate and other factors (Campbell et al., 2017; Hou et al., 2022; Lavergne et al., 2022; Zscheischler et al., 2014) are at least as important as vegetation area for GPP simulation. Therefore, achieving consistency in the area for vegetation types is not sufficient, and an improved understanding and representation of C₃ and C₄ GPP in models are needed.

The range in total area fraction of C₄ vegetation of 10 %–26 % in the CMIP6 models spans the observation-based estimates of 17.1 %–17.7 % (Still et al., 2003; Luo et al., 2024) and has a similar range as the dynamic global vegetation models in the TRENDY project (7 %–23 %; Luo et al., 2024). Since the observation-based estimates are quite consistent, model simulations of carbon fluxes can likely be improved by detailed comparison and specification of C₄ vegetation cover in particular PFTs within models to better align them with observational estimates.

The trends in C₃ vs. C₄ contributions to GPP simulated by CMIP6 models produced differing trends in stable carbon isotopic discrimination when only C₃ and C₄ fractions of GPP are used to calculate discrimination trends Eq. (1) (Fig. 6). The simulated changes in C₃ and C₄ fractions of GPP did not produce strong increasing trends as found by analysis of atmospheric δ¹³C data (Keeling et al., 2017; Peters et al., 2018). Since the effect of changing C₃ and C₄ fractions of GPP on stable carbon isotopic discrimination varied widely across the CMIP6 models, it is still uncertain how much they could add to or oppose changes in stable carbon isotopic discrimination caused by environmental or physiological effects. Nevertheless, Lavergne et al. (2026) found that accounting for the changes in C₃ and C₄ vegetation abundance improves the modelled trends in atmospheric carbon isotope composition compared to the data. In particular, since the GPP trend in the models (11 % to 18 %, Fig. S4) is weaker than the CO₂ fertilisation effect over the 20th century based on carbonyl sulphide data (+ 30 %; Campbell et al., 2017), there may have been a more positive C₃ vegetation-driven trend in stable carbon isotopic discrimination in reality than in the models. Since the CMIP6 models did not provide output for GPP for C₃ and C₄ vegetation or for crops and natural vegetation, we had to scale total GPP by the area fractions of each vegetation type for grid cells with mixed vegetation. This calculation may not represent the models' actual attribution precisely, since C₄ vegetation can have higher GPP relative to area (Luo et al., 2024) in grid cells with a mosaic of C₃ and C₄. However, the extra resources needed in the C₄ photosynthetic pathway (Ehleringer, 1978) to concentrate CO₂ may make the difference between C₃ and C₄ CO₂ uptake minimal, especially if they co-occur relatively equally in such grid cells. Incorporating more detailed vegetation information from the models would improve the accuracy of the calculations, and analysis based on TRENDY models could be a good place to start.

5 Conclusions

We analysed changes in vegetation with C₃ and C₄ photosynthetic pathways for natural vegetation and crops in 11 CMIP6 models over the historical period of 1850 to 2014. Except for one model that used a fixed vegetation distribution, the models include the expansion of agriculture using LUH2 data. Still, there is a significant variation in the fraction of area and GPP allocated to crops in these models, and the UKESM1 model strongly underestimated crops in Asia and Africa, likely due to biases in simulated climate. The C₄ fraction of vegetated area has remained relatively constant in the models, though there is a large range in simulated C₄ area fraction of 10 % to 26 % in 2014, while observation-based estimates are quite consistent at 17.1 %–17.7 % from the 1980s to 1990s and between 2001 and 2019 (Luo et al., 2024; Still et al., 2003). Overall, the total vegetation area in most models is not changing, but the global vegetation composition is changing in favour of crops. Whereas all the models agree that C₃ vegetation GPP and total global GPP are increasing, the magnitude of the increase spans a wide range. In C₄ vegetation, the models disagree on the magnitude of GPP and the sign of its trend, especially in natural vegetation, implying that our understanding of and model parametrisations for the simulation of C₄ GPP need improvement. The strong positive trends in UKESM1 C₄ crop and MPI-ESM-1-2-HAM and MPI-ESM1-2-LR C₄ natural vegetation area clearly drive the decline in their global stable carbon isotopic discrimination trend, an effect that is not as pronounced in other models' vegetation categories. However, due to the large uncertainty in vegetation area abundance and GPP, it will be difficult to determine if this effect of change in vegetation area abundance and GPP on the global stable carbon isotopic discrimination trend is consistent across all the models. We have shown in this study that carbon flux simulations in current ESMs include uncertainties in the C₃ and C₄ photosynthetic pathways in crops and natural vegetation; reducing these uncertainties will require stronger observational constraints and more robust and realistic model parameterisations.