Modelling approach
We used the optimality-based light-use efficiency P model19–21 to simulate optimal ecosystem gross primary production and carbon isotopic discrimination and develop a simple approach based on the P model simulations to predict the share of C4 plants in the total GPP (see Text S1 for more details on the model and workflow in Figure S1). We converted the potential share of C4 plants in the total GPP to fraction of C4 plants (F4,pot) using the emergent constraint from Luo et al. (2024)6 as the share of C4 plants divided by 1.13.
Since the total fractions of C3 and C4 plants consist of natural ecosystems (trees, grasses) and crops, we estimated the C4 fraction of natural ecosystems by removing the fraction of human managed areas (C3 and C4 crops and urban areas) from F4,pot as:
$${F}_{4, natural} = {F}_{4,pot}-\left({F}_{3, crops}+{F}_{4, crops}\right){-F}_{urban}$$
1
We then calculated the total fraction of C4 and C3 plants (F4,tot and F3,tot) considering natural ecosystems (F4,nat and F3,nat) and crops (F4,crops and F3,crops) as:
$${F}_{4, tot} = {F}_{4,nat}+{F}_{4, crops}$$
2a
$${F}_{3, tot} = {1-F}_{4,tot}= {F}_{3,nat}+{F}_{3, crops} \left(2b\right)$$
We estimated C3, C4 and total (C3 + C4) GPP from their respective potential GPP as:
$${GPP}_{C3}={GPP}_{C3,pot}{F}_{3,tot}$$
3a
$${GPP}_{C4}={GPP}_{C4,pot}{F}_{4,tot}$$
3b
$${GPP}_{tot}={GPP}_{C3}+{GPP}_{C4}$$
3c
We aggregated GPP for C3, C4 and all (C3 + C4) plants from each grid cell weighted by grid-cell area for each prediction with each C4 map and compared their respective temporal changes (GPP, in PgC yr-1).
We then estimated the stable carbon isotopic discrimination (Δ) as:
$${\Delta }={F}_{4,tot} {{\Delta }}_{4}+{F}_{3,tot}{{\Delta }}_{3}={{\Delta }}_{\text{n}\text{a}\text{t}\text{u}\text{r}\text{a}\text{l}}+{{\Delta }}_{\text{c}\text{r}\text{o}\text{p}\text{s}}$$
4
with Δ4 and Δ3 the carbon isotope discrimination of C4 and C3 plants, respectively (see Text S1), calculated for each month (t) and weighted by each month’s C3 and C4 GPP to obtain annual averages of Δ3 and Δ4:
$${\varDelta }_{yr,x} =\sum _{t}\frac{{\varDelta }_{x,t}\times {GPP}_{x,t}}{\sum _{t}{GPP}_{x,t}}$$
5
with x being the plant type (C3 or C4) selected.
Model configuration and simulation
We run the P-model on a 0.5° resolution grid at a monthly timestep over the period 1982–2016, forced by monthly mean values of daytime air temperature (Tair, ºC) and vapour pressure deficit (VPD, Pa), incident photosynthetic photon flux density (PPFD, mol m–2 month–1), the fraction of incident PPFD absorbed by foliage (fAPAR, dimensionless), the atmospheric partial pressure of CO2 (ca, Pa), root-zone volumetric soil moisture (\({\theta }\), m3 m–3) and elevation (z, m). Monthly mean Tdaytime and VPD were calculated from the Climatic Research Unit gridded time-series (CRU TS4.03) dataset37 to consider only the part of the day when photosynthesis occurs, as described in Lavergne et al. (2020)38. Monthly shortwave downwelling radiation (SWdown) was obtained from the daily WATCH-Forcing-Data-ERA-Interim (WFDEI) data39 and used to calculate monthly PPFD.
Monthly fAPAR data were derived from the Advanced Very High Resolution Radiometer (AVHRR) Global Inventory Modeling and Mapping Studies (GIMMS) fAPAR 3g product40 gridded at 0.5° resolution. Since monthly ca can vary strongly across the world, we used the annual ca data (µmol mol–1) derived from 41 and converted them into Pascal using z derived from WATCH-WFDEI. Monthly \({\theta }\) (m3 m–3) over a 1 m soil depth was calculated using a modified version of the SPLASH model42 driven by daily precipitation, Tair and SWdown from the WATCH-WFDEI dataset over 1979–2018. The model \({\theta }\) outputs were derived from 43.
Annual percentage of treecover used to estimate the fraction of C4 plants (see Text S1) were derived from NASA Making Earth System Data Records for Use in Research Environments (MEaSURES) Vegetation Continuous Fields (VCF)5KYR v00144 for the 1982–2016 period. The years 1994 and 2000 were missing from this dataset, we therefore interpolated treecover for these years averaging values from the previous and subsequent years (1993 and 1995, and 1999 and 2001, respectively).
We used the urban areas and C3 and C4 crop distribution estimates from the LUHv2-2019 dataset (https://daac.ornl.gov/VEGETATION/guides/LUH2_GCB2019.html). Data for the remote-sensing products given at 0.05° resolution were aggregated to 0.5° resolution using the mean of all the 0.05° grid cells within the 0.5° grid cell.
Model evaluation and comparison
We evaluated the simulations of Δ13C for C3 and C4 plants using a compilation of leaf isotopic data22 that includes 3601 measurements for C3 plants and 531 measurements for C4 plants (see Figure S4a for site locations). We also evaluated the model predictions of the fraction of C4 plants (F4) over the 1982–2016 period using a new compilation of 2156 soil δ13C measurements derived from published sources specifically compiled for this study23 (see Table S1 and Figure S4b for data information and locations). Predicted Δ13C values for C3 and C4 plants were converted to δ13C using atmospheric δ13CO2 from Graven et al. (2017)45 to enable comparisons with the observed network.
We compared the predictive skills of our approach with those from two independent, global C4 distribution maps. We used the widely-used map from Still et al. (2009)14 available at https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=932 and the map recently developed by Luo et al. (2024)6 for the period 2001–2019 (https://zenodo.org/records/10516423) - referred to as Still2009 and Luo2024, respectively. Since the studied maps covered different regions of the world, we homogenized them to make them comparable. For instance, when grid points were set to "NA" but were in land surface areas (excluding deserts and snow/glacial lands), we assumed F4 was 0. As a result, the estimated F4 values derived from the maps may be slightly different from those highlighted in the respective publications. We predicted soil isotopic composition using each of the maps and compare them with the soil isotopic network. We also compared the mean averages and trends in total F4, GPP and Δ13C derived from each of the maps.
All model runs, statistical analyses and figures were performed using Python.