Vegetation biogeography is a main source of uncertainty in modelling the land carbon cycle

Domain: biology
Published: 2025-12-17
Editorial outcome: minor_revision

Abstract
The terrestrial biosphere exchanges a large amount of CO2 with the atmosphere through photosynthesis and respiration, determining the magnitude of land carbon sink and consequently influencing the rate of global warming. The magnitudes of global photosynthesis and respiration, however, vary widely across models (100-200 PgC/year), constituting a key and persistent source of uncertainty in carbon cycle and climate modelling. Here, we argue that the uncertainty in the land carbon cycle modelling is largely attributable to the uncertainty in biogeography – the distribution of plant functional types (PFTs). Using an ensemble of dynamic global vegetation models (DGVMs), we find a strong dependence of total photosynthesis on total area for each PFT. The dependence allows us to reduce the spread of land carbon cycle estimates by ~75% using remote sensing-based PFT maps. We further find that 56 ± 21% of climate-driven changes in global photosynthesis modelled by DGVMs are caused by changes in PFT distribution in the last two decades. Our study identifies vegetation biogeography as a main controlling factor of uncertainty in land carbon cycle modelling and highlights the importance of biogeography-climate interactions in carbon cycle and climate studies.

Introduction
The land carbon cycle is a critical component of the global carbon budget1, modulating the warming rate of the global climate by absorbing approximately 25% of anthropogenic CO2 emissions2,3. Photosynthesis (i.e., gross primary productivity; GPP) and ecosystem respiration (ER) are the two primary carbon fluxes in the land carbon cycle and determine the magnitude of the sink1. However, global estimates of GPP and ER from models are highly uncertain, ranging from roughly 100 to 200 PgC/year4,5. The considerable uncertainty in the land carbon cycle persists across several generations of dynamic global vegetation models (DGVMs) and Earth system models (ESMs), posing a long-standing challenge in carbon cycle simulations and climate predictions6,7.

Considerable efforts have been made to identify the sources of uncertainty in models for land carbon cycle simulations, such as assessing climate forcing datasets8–10 and improving representations of ecological processes (e.g., thermal acclimation of photosynthesis11, explicit representation of microbial processes12, plant hydraulics13, and competition for resources14), and parameterizations (e.g., photosynthetic capacity15,16 and microbial carbon use efficiency17) in process-based models. However, these adjustments have yet to significantly reduce the uncertainty in GPP and ER, as evidenced by the similar spread of land carbon simulations from DGVMs over the past decades1,18. Meanwhile, remote sensing-based proxies19 and data-driven machine learning approaches20,21 have also been developed to provide estimates of GPP and ER. While these approaches are reported to estimate global GPP and ER from observations and have less uncertainty, they have limitations in identifying the mechanisms that are missed or misrepresented in carbon cycle simulations. These limitations include prescribed climate sensitivity22, lack of representation in the tropics20, empirical nature of the scaling factor between GPP and spectral signals23, and lack of or inconsistent consideration of CO2 fertilization effect24,25.

Here, we hypothesize that the key source of uncertainty in land carbon cycle simulations is biogeography–the distribution of plant functional types (PFTs), based on multiple lines of evidence. Biogeography is known to affect the simulations of the land carbon cycle in DGVMs26–28, causing up to a 24% difference in global GPP simulations within a single model26–28. The large impact implies that the inter-model differences in estimated land carbon budget are potentially linked to the inter-model differences in biogeography simulations. In particular, a recent study29 shows that the inter-model differences in C4 vegetation distribution, a key PFT, explained over 80% of the uncertainties in GPP estimates for C4 vegetation, providing direct evidence of the role of biogeography in land carbon cycle simulations. Additionally, the current framework for improving land carbon cycle simulations in DGVMs implies the importance of biogeography, given that changes in model parameters, structures, and even the emerging climate sensitivities of simulated carbon fluxes can also be driven by PFT change, as DGVMs largely regard these processes as PFT-specific26,30–32.

To test the hypothesis on the dependence of the simulated land carbon cycle on biogeography, we use the outputs from an ensemble of state-of-the-art DGVMs to assess the relationships between total GPP and total area for each PFT. We aim to quantify the dependence of PFT total GPP on PFT distribution and quantify the contribution of vegetation biogeography to the uncertainties in global GPP and ER estimation (see Methods). The dependence will then be applied to narrow inter-model discrepancies in GPP and ER using remote sensing-based (RS) PFT maps. Using the PFT-specific relationships between total GPP and total area, and the simulated changes in PFT distributions from some DGVMs, we further explore the contribution of biogeography changes to the changes in global GPP over the past two decades.

Results
The relationships between PFT GPP and area

We first established the relationships between total GPP (i.e., in PgC/year) and total area (i.e., in km2) for each PFT, using the simulations from an ensemble of DGVMs (see Methods). For each of the twelve PFT types on the land surface, we observed a strong and significant relationship (p < 0.001) between GPP and area across models (Fig. 1). Using alternative approaches (i.e., linear regression without forcing through zero or multivariate regression), we consistently observed the strong dependence of total GPP on area for most of the PFTs (Fig. 1; Supplementary Fig. 10; Supplementary Table 5). The strong relationships indicate that the inter-model discrepancy of total GPP of each PFT is highly relevant to the difference in the prescribed or simulated PFT area across DGVMs, and that DGVMs simulate a similar per-area GPP on average for each PFT (i.e., the slopes of the GPP-area relationships in Fig. 1). Specifically, the average per-area GPP is 2.48 ± 0.07 PgC/year per 10^6 km^2 for evergreen broadleaf forests (EBF), 1.38 ± 0.08 PgC/year for temperate evergreen needleleaf forests (ENFt), 0.94 ± 0.02 PgC/year for boreal evergreen needleleaf forests (ENFb), 1.22 ± 0.01 PgC/year for deciduous broadleaf forests (DBF), 1.10 ± 0.02 PgC/year for deciduous needleleaf forests (DNF), 0.66 ± 0.02 PgC/year for shrubs (SHR), 0.78 ± 0.06 PgC/year for warm C3 grasses (C3Gw), 0.67 ± 0.03 PgC/year for cool C3 grasses (C3Gk), 0.34 ± 0.04 PgC/year for cold C3 grasses (C3Gc), 0.94 ± 0.10 PgC/year for C3 crops (C3C), 1.40 ± 0.09 PgC/year for C4 grasses (C4G), and 1.38 ± 0.04 PgC/year for C4 crops (C4C).

Updated land carbon cycle estimation using remotely sensed PFT area

Since there is a large spread in PFT area across DGVMs, the uncertainty in PFT area propagates to the substantial uncertainty in global GPP estimation (Fig. 1). The global GPP estimated by the ensemble of DGVMs ranged from 113.7 to 194.0 PgC/year, with a mean and standard deviation (SD) of 146.32 ± 26.39 PgC/year (Fig. 2a). The global ER estimated by the DGVMs was 143.39 ± 25.99 PgC/year (Fig. 2b), and had the magnitude similar to GPP since global GPP and ER were strongly coupled across DGVMs (Fig. 2c). The tight relationships between PFT-specific total GPP and area across DGVMs (Fig. 1) allowed us to reduce the uncertainty in global GPP and then ER (Fig. 2c) if we have robust PFT distribution maps from observations.

We, therefore, used the PFT area from five widely-used RS PFT maps – including MODIS v6.133, ESA CCI v2.0.834, CGLS_LC100 v335, GLC_FCS30D36, and GLASS-GLC37 (see Methods) – in the GPP-area relationships (Fig. 1) to estimate the PFT-specific total GPP, and then summed up to global GPP. We found that the global GPP was updated to 135.66 ± 6.54 PgC/year on average from 2001 to 2019. The uncertainty in the global GPP estimate was reduced by approximately 75% (i.e., from 26.39 PgC/year to 6.54 PgC/year; Fig. 2a). The reduced uncertainty in GPP is mainly due to the reduced uncertainty in RS PFT mapping compared to those used in DGVMs. Specifically, the total GPP of PFTs are 41.35 ± 4.67 PgC/year for EBF, 6.76 ± 1.41 PgC/year for ENFt, 5.20 ± 0.92 PgC/year for ENFb, 14.95 ± 2.59 PgC/year for DBF, 5.27 ± 0.93 PgC/year for DNF, 8.21 ± 1.84 PgC/year for SHR, 7.51 ± 1.03 PgC/year for C3Gw, 5.80 ± 1.08 PgC/year for C3Gk, 1.42 ± 0.20 PgC/year for C3Gc, 10.55 ± 2.36 PgC/year for C3C, 24.07 ± 0.98 PgC/year for C4G, and 4.51 ± 0.15 PgC/year for C4C. The constrained global GPP was within the range of other GPP estimates such as those based on solar-induced fluorescence19,38 and soil respiration39, but lower than those inferred from plant carbonyl sulfide40 and higher than those upscaled from eddy covariance observations using machine learning20,21 (Fig. 2a). The constrained global GPP was further used to estimate global ER (Fig. 2c), with the uncertainty reduced from 25.99 to 6.53 PgC/year (Fig. 2b). Using linear regression not forced through zero, we still achieved a substantial yet smaller reduction (i.e., 48% compared to 75%) in the spread of inter-model GPP and ER (Supplementary Fig. 11).

Contribution of PFT change to global GPP change

We further used the relationships between PFT-specific total GPP and area to estimate the changes in global GPP induced by PFT area changes (Fig. 2d), using the annually adjusted PFT area from four DGVMs that enabled dynamic biogeography simulations (Supplementary Table 1). From 2001 to 2019, these DGVMs demonstrated that elevated CO2 increased global GPP by 23.47 ± 3.25 PgC/year and climate change increased global GPP by 2.05 ± 0.45 PgC/year. Of the global GPP increase induced by elevated CO2, we found that approximately 20 ± 4% of the increase resulted from elevated CO2-driven PFT changes, primarily due to a decrease in C4 grass distribution and an increase in C3 plant distribution simulated by these DGVMs (Supplementary Fig. 3). Similarly, of the global GPP increase induced by climate change, we found that 56 ± 21% resulted from climate change-driven PFT changes. Among these DGVMs, the climate change-driven PFT changes include the expansion of C4 grasses, C4 crops, and deciduous forests in high latitudes, in replacement of C3 grasses and evergreen forests. The expansion of C4 was likely due to its anatomical traits, which confer higher water-use efficiency41,42.

Large uncertainty in biogeography simulations

We further examined the multi-year average PFT distribution simulated by each DGVM in the framework of Whittaker bioclimatic schemes43 (Supplementary Fig. 4), where the distribution was presented in a space of mean annual temperature (MAT) and mean annual precipitation (MAP). While these biogeography simulations showed similar PFT types along the gradient of warm to cold (e.g., EBF to DBF, and to ENF), and wet to dry (e.g., EBF to SHR, and to C3G and C4G), there were clear differences in the bioclimatic boundaries of all PFTs between models, as well as between models and RS PFT maps. For example, the MAT range for C4G was narrow (20–30°C) for CLASSIC and RS PFT maps, but much broader (-20–30°C) for JULES. The MAP range for DBF was 0–2000 mm/year for JSBACH and 1000–2000 mm/year for RS PFT maps, but could reach 0–4000 mm/year for LPX_Bern and ORCHIDEE (Supplementary Fig. 4). Examining the spatial distribution of PFTs across models (Supplementary Fig. 5; Supplementary Fig. 6), we identified the hotspots for PFT disagreement between DGVMs. Those hotspots of biogeography uncertainty include EBF in the southern U.S. and Mexico, ENF in the boreal-temperate transition zones, DNF in North America, C3G in the low latitudes, C4G in Central Asia, and DBF and SHR across the globe. Given the difference in the mechanisms each model uses to adjust PFT distribution, it is challenging to generalise the reasons for biases in PFT distribution across models.

Moreover, we found that RS PFT maps also show large differences in their estimates of PFT total area and spatial distributions (Fig. 3; Supplementary Table 3). For example, there were large spreads in the estimates of global EBF area (i.e., from 13.97 to 21.56 × 10^6 km^2) and C3C area (from 7.66 to 16.34 × 10^6 km^2) across the five remote sensing products. The hotspots of biogeography disagreement among RS PFT maps include EBF in the eastern U.S. and western Europe, ENF in northern Canada and eastern Russia, DNF in North America and Europe, and SHR and C3G across the globe. The discrepancy across the RS PFT maps propagated into the uncertainty of our constrained global GPP (Fig. 1; Fig. 2).

Methods
Dynamic global vegetation models (DGVMs)

To examine the relationship between the PFT total GPP and the corresponding PFT area, we used simulations from an ensemble of DGVMs participating in the TRENDY v9 project72. The DGVMs were coordinated to perform simulations under four different scenarios (S0, S1, S2, and S3), following the TRENDY protocol (https://blogs.exeter.ac.uk/trendy/protocol/). S0 is the control simulation, with no changes in the external forcings for DGVMs (i.e., time-invariant pre-industrial CO2, climate and land use mask). In S1, only atmospheric CO2 concentration varied when forcing DGVMs. In S2, atmospheric CO2 and climate varied. In S3, atmospheric CO2, climate and land use varied when forcing DGVMs.

We mainly used the simulations from the S2 scenario to examine the relationship between PFT-specific total GPP and area. Changing the simulations to other scenarios, such as S3, would not affect the relationship between PFT total GPP and PFT area across models (Supplementary Fig. 9), however, we would like to exclude the consideration of land use and land cover change in this analysis. A total of 11 models with varying spatial resolutions provide the required PFT-specific outputs (Supplementary Table 1). We obtained the annual total GPP and total area for each PFT for the period 2001 to 2019, and examined the relationships between them using the multi-year averages. Since not all models adopted the same PFT classification system, we harmonized the PFT classifications in each DGVM into twelve common PFT types, which cover the major functional differences in plant photosynthetic pathways, climate adaptation, and phenological cycles. These PFTs include evergreen broadleaf forests (EBF), temperate evergreen needleleaf forests (ENFt) and boreal evergreen needleleaf forests (ENFb), deciduous broadleaf forests (DBF), deciduous needleleaf forests (DNF), shrubs (SHR), warm C3 grasses (C3Gw), cool C3 grasses (C3Gk), and cold C3 grasses (C3Gc), C4 grasses (C4G), C3 crops (C3C), and C4 crops (C4C). The availability of PFT types in each DGVM is presented in Supplementary Table 1. For those DGVMs that only provide evergreen needleleaf forests (ENF) and C3 grasses (C3G) without subtypes, we used a set of bioclimatic boundaries73 to divide ENF into ENFt and ENFb, and to divide C3G into C3Gw, C3Gk, and C3Gc (Supplementary Table 2), as distinct climates otherwise might mask the relationship between PFT area and GPP (Supplementary Fig. 1). The climate data is obtained from CRUJRA v2.174, which is also the climate forcing data used for DGVMs in the TRENDY v9 project.

We also used PFT-specific GPP outputs from DGVM, which was used to establish PFT-specific relationship between total GPP and total area (Fig. 1). However, since DGVMs under TRENDY v9 did not supply PFT-specific ER, we were unable to establish PFT-specific relationship between total ER and area. We were only able to establish a relationship between global GPP and global ER in our analysis (Fig. 2c).

Remote sensing estimates and uncertainties of global PFT distribution

We adopted five widely-used remote sensing-based (RS) PFT products, including the MODIS33, ESA CCI34, CGLS_LC10035, GLC_FCS30D36, and GLASS-GLC37, to constrain the GPP estimates based on the relationship between PFT-specific total GPP and total area.

The MODIS and ESA CCI datasets directly provide PFT maps (MCD12Q1v6.1 and ESACCIv2.0.8), with some assessments of their production uncertainties28,65. Specifically, the MODIS PFT map was translated from land cover classification within the same product at a spatial resolution of 500 metres using a series of climate-based rules33,75. ENFt and ENFb were classified from the ENF land cover type, with ENFt requiring a minimum temperature of the coldest month (Tc) above -19°C and an accumulated growing degree-days above 5°C (GDD) exceeding 1200°C, while ENFb had colder conditions75. C3 and C4 grass were distinguished within grassland cover with GDD > 1000°C, Tc > 22°C and monthly precipitation > 25 mm as C4 grass, others as C3 grass75. C3 grass was further classified to C3Gw, C3Gk, and C3Gc subtypes using DGVM rules (Supplementary Table 2). EBF, DBF, DNF, and SHR were mapped from their corresponding land cover types from MODIS33.

The ESACCI PFT map was translated from the ESACCI land cover product at a 300-metre resolution from 1992 to 2020, with high-resolution auxiliary datasets such as tree cover and tree canopy height data34,76,77. Pixels classified as tree types (e.g., EBF, DBF, ENF, and DNF) in the land cover maps were assigned to corresponding PFTs, with tree cover percentage determined using a 30-metre resolution tree cover dataset. For shrubland classes (a mix of tree and shrub woody vegetation types), a tree canopy height dataset was used to separate shrubs from trees34,75.

The CGLS_LC100 provided the annual fraction of land cover classes, including EBF, ENF, DNF, DBF, SHR, mixed forest (MF), unknown forest (UNF), cropland and grassland, at a 100-m resolution from 2015 to 201935. The GLC_FCS30D provided the categorical information for the same classes at 30-m grid from 1985 to 2022, updated every five years before 2000 and annually thereafter36. The GLASS-GLC provided annual distribution of forest (F), SHR, grassland and cropland at 5-km resolution from 1982 to 201537. After aggregating to 0.5-degree resolution for each land cover map, the MF fractions were translated to forest PFTs and SHR according to the translation rules (Supplementary Table 2). The UNF and F were translated to EBF, ENF, DBF, and DNF PFTs according to the same criteria. To make the PFTs map compatible with those from DGVMs, we further divided the RS PFTs using environment constraints78 when needed (i.e., such as from ENF to ENFt and ENFb).

Since all the RS PFT maps lack explicit information on C3/C4 grass and C3/C4 crop classification, we used a newly developed C4 vegetation map30 to divide croplands and grasslands into their C3 and C4 counterparts. The C4 vegetation map was generated using a global database of photosynthetic pathways, satellite observations, and a photosynthetic optimality theory. It contains C4 grass and C4 crop distribution over the globe at a 0.5-degree resolution from 2001 to 201929.

We quantified the uncertainty in RS PFT maps as the one standard deviation of PFT area among MODIS, ESA CCI, CGLS_LC100, GLC_FCS30D and GLASS-GLC maps. We did not consider the production uncertainty of individual RS PFT maps as there was a lack of assessments for most products (Supplementary Table 3)28,65.

The relationship between PFT total GPP and PFT area

In this study, we identified the relationships between the simulated total GPP and total area for each PFT across DGVMs (i.e., GPP-area relationships). We obtained the multi-year average total GPP and total area for each PFT over the period 2001-2019, and examined the relationships between them.

GPP_PFT_i = β_i × area_PFT_i (1)

where, β_i is the relationship between total GPP (GPP_PFT_i) and area (area_PFT_i) across DGVMs for i^th PFT. The i ranged between 1 to 12, representing EBF, ENFt, ENFb, DNF, DBF, SHR, C3Gw, C3Gk, C3Gc, C4G, C3C, and C4C, respectively. We forced the intercepts of the linear models to zero, based on the assumption that PFT-specific GPP should be zero if the PFT area is zero. Meanwhile, we also examined the relationships without this assumption. We estimated the uncertainties in β_i through a bootstrapping method73,79, where one of the DGVMs was removed each time for regression to get the upper and lower boundary of the relationships.

Based on the GPP-area relationships and area of each PFT from RS PFT maps (RSarea_PFT_i), we estimate the constrained global GPP (GPP_global_C) over the globe.

GPP_global_C = Σ_{i=1}^{12} β_i × RSarea_PFT_i (2)

We also identified a relationship between global ecosystem respiration (ER) and global GPP across DGVMs,

ER_global = γ × GPP_global (3)

where, γ is the slope, ER_global and GPP_global stand for the global ER and global GPP in each DGVM. The global ER in each model is calculated by the sum of the available autotrophic respiration and heterotrophic respiration. Combining γ and the constrained global GPP (GPP_global_C), we estimated the constrained global ER (ER_global_C).

ER_global_C = γ × GPP_global_C (4)

Machine learning emulators

To assess the relative importance of biogeography and biogeochemistry in GPP modelling, we developed a series of emulators for DGVMs by training random forest models globally for each PFT of each DGVM. Specifically, each emulator was trained to predict PFT GPP rate on grid cell level using annual climate forcings (temperature, precipitation, radiation, pressure, and wind speed) from the CRUJRA v2.1 datasets74 from 2001 to 2019, which was consistent with the climate forcing data for DGVMs in the TRENDY v9 project. We excluded the LPJ-GUESS model from our emulator analysis, as in this model, the pixel PFT GPP rate is not available to support the training. We conducted the training process in MATLAB R2024b, with automatic tuning of key parameters, including the number of trees and the minimum leaf size. These emulators were independently validated via cross-validation, demonstrating consistently high predictive accuracy (mean R² values > 0.8 for all PFTs, Supplementary Table 4). After developing these emulators, we used them to estimate PFT-specific GPP rate within each grid cell and obtain PFT-specific GPP by using a consistent PFT fraction map from MODIS33. We then aggregated the PFT-specific GPP of all grid cells and calculated global GPP.

The contribution of PFT changes to GPP changes

Using the relationships between PFT-specific GPP and area, we also estimated the contribution of elevated CO2- and climate change-driven PFT changes to the GPP changes. Specifically, we first obtained the simulated PFT area from DGVMs under S0 (time-invariant CO2 and climate), S1 (varying CO2 and time-invariant climate), and S2 (varying CO2 and climate) scenarios. By subtracting the PFT areas of S0 from S1, we calculated the PFT area changes driven by elevated CO2. Similarly, by subtracting the PFT areas of S1 from S2, we achieved the changes in PFT driven by climate change. Then we used the PFT area changes in the GPP-area relationship (Fig. 1) to infer the associated change in GPP. These PFT-driven changes in GPP were further compared against the total GPP changes driven by elevated CO2 and climate change. Please note that in this analysis, we only used the four DGVMs (i.e., JULES, LPJ-GUESS, LPX-Bern, and SDGVM) that dynamically adjust PFT distribution over time (Supplementary Table 1).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Land carbon sink remain largely uncertain among intercomparison models, thus restricting our capability of predicting, mitigating and adapting to climate change. There have been research efforts which leveraged the emergent constraint approach to reduce the predicted uncertainty for the variables (i.e., photosynthesis or gross/net primary production - GPP/NPP and ecosystem respiration – ER). The essence of emergent constraint is to identify the emergent relationships (typically linear) between a multi-model ensemble and an observational estimate. Here this study by Dr. Zhao et al have used this emergent constraint approach to identify the remarkably strong linear relationship between carbon flux (GPP and ER) and biogeographic area of vegetation biomes. Then the two sets of remote sensing datasets (MODIS and ESA CCI PFT maps) were used as the observational estimate to constraint the uncertainty of GPP and ER, respectively. I am highly impressive by the approach of using biogeographic vegetation biome area and its strong relationships because this offers significant insights in term of vegetation biome impact on carbon cycling.\nHere I have some moderate suggestions to improve this study. First, the key uncertainty or weakness (if any) is that the remote sensing data of vegetation biome area has large uncertainty as shown in the results of figure 1 et al. The authors have made a good job to clarify this. Since the emergent constraint largely depends on the accuracy of observational variable – vegetation biome area here. The authors may consider clarifying that 1) the emergent constraint could still reduce the uncertainty even if the two sets of remote sensing datasets were used together by combing the remote sensing datasets (i.e., by getting the range of minimum and maximum values of GPP and ER); 2) the needs for future studies to improve or have more reliable observational estimate of vegetation biome area. Second, the authors may also clarify the potential limitation of emergent constraint in term of underlying mechanisms underlying mechanisms (i.e., Winkler et al., 2019) and the spatial/temporal scales mismatch between data and models (i.e., Winkler et al., 2019). Ideally or necessarily, the research efforts in reducing the uncertainty of land carbon sink should be aimed by both improving the model structure or processes (i.e., the emergently vegetation biome distribution) and the emergent constraint approach here. This also related to the constraint of ER in this study because the current model structure have very simple structure for modeling ER which are strongly coupled with GPP. Third, the authors may also consider the potential effects of land usage impacts. Land use change likely have influenced the vegetation biome distribution more than climate change in term of spatial and temporal scales. Currently, it is not clear how the land usage change has been accounted for.\nSince the paper is very well written with insightful discussion, I also have the following minor suggestion for the authors to consider\nAbstract: may also consider reporting the results of ER\nLines 44-45: it would be important to mention competition, i.e., by accounting for size or cohort based competition for resources (light, nutrients and water)\nIntroduction: may consider clarifying what the emergent constraint is and the limitation of emergent constraint.\nFigure 1: the color scheme would need to be improved. Also there are different definitions for vegetation biomes such as WWF vegetation biome. My knowledge is also that there is 1 km resolution of land cover maps https://www.earthenv.org/landcover I am not sure whether this is useful as an additional observational constraint.\nOverall, the whole manuscript has focused on the constraint of GPP, but the results of ER was indeed also important for carbon sink. Thus, the authors may consider balancing this.\nFigure 2: I am curious whether N deposition has been accounted for and how important it would be.\nDiscussion: as clarified, the authors may consider clarifying the limitation of emergent constraint and needs of improving model structure/processes for the mechanistically deriving/model the emergent properties of vegetation biome.\nThe section of discussion on structure (vegetation biome biogeographical distribution) vs biogeochemistry is interesting. It is also good to learn that the authors used emulator to disentangle their relative impacts. With this, the authors may consider elaborating a bit more for the implication of improving the model mechanisms in term of more research efforts on vegetation biome biogeographical distribution itself vs biogeochemistry.\nReferences\nWinkler, A. J., Myneni, R. B. & Brovkin, V. Investigating the applicability of emergent constraints. Earth Syst. Dyn. 10, 501–523 (2019).", "reviewer_2": "Reviewer’s comments\nThis paper attempts to identify an observational constraint for dynamic global vegetation models (DGVMs) by establishing a linear regression relationship between plant function type (PFT) area and total gross primary productivity (GPP) for a given PFT. The authors use satellite-based observations to infer a plausible range of total global GPP per-PFT (GPPpft). They correctly identified that land surface models with dynamic vegetation components are often deficient when it comes to simulating present-day biogeography, and that using these PFT observations to constrain total global GPP estimates can tighten the range of GPP uncertainty.\nI have serious concerns, however that the linear regression on which this study is based is fundamentally flawed because both variables used in the regression have a strong dependence on PFT area. While I broadly agree with the conclusion that observed global area of a PFT can be used to refine global GPP estimates derived from DGVMs, I have serious reservations about the approach taken here. My criticism boils down to the fact that in this regression relationship, the independent variable (x-axis; PFT area) is also a central component of the dependent variable (y-axis; total global GPPpft), and therefore the resulting strong relationships shown for each PFT are not valid because it is essentially using PFT area to predict a variable that is strongly dependent on PFT area. In fact, the only reason why there is variation in the relationship (and not a perfect 1:1 line) is because there is variation in the per-model and per-grid cell GPPpft rate per unit area. In order to derive total global GPPpft, I assume the authors would have had to multiply the raw land surface model output GPPpft rate per unit area (likely in units similar to gC m-2 yr-1) by the grid cell PFT area, and then sum all grid cells (possibly also making a multi-year mean of the timeseries). Similarly, converting model output of PFT fraction to a total land area for a given PFT (Areapft) would involve multiplying the per grid cell fraction by each grid cell’s area (giving the exact value used in the GPP calculation in the previous sentence), and then summing all grid cells. If you take out the area dependence from the simple GPP calculation, you are left with the GPPpft rate again per unit area. Therefore, I think the linear regression on which the emergent constraint is based is flawed.\nHaving said that, I do think there is value in highlighting that observations of PFT area can be used to constrain the simulated GPP of models in order to gain a more accurate measure of present-day GPP and consequently atmosphere-land carbon fluxes. I just don’t think that it can be done in the way that the authors have proposed. I believe that the authors would learn more about the models and the plausible constraints on GPP by looking in detail at the locations where PFT fractions match observations, and those where PFT fractions do not. I would also encourage the authors to examine other carbon-cycle emergent constraint studies and consider the definition presented in the Cox et al. (2019) paper:\n“The emergent constraint technique involves using the full ensemble of models to find an across-ensemble relationship between an observable feature of the Earth System (such as a trend, interannual variation or change in seasonality) and an uncertain aspect of the future.”\nI would argue that while PFT fraction is an observable feature of the Earth System (and GPP the uncertain aspect of the future), I don’t think that it provides a trend, interannual variation or change in seasonality, in the same way that other emergent constraint studies have. See for example, Cox et al. (2013) who found a strong linear relationship between the sensitivity of tropical land carbon storage to warming and the sensitivity of the annual growth rate of atmospheric CO2 to tropical temperature anomalies.\nOther Comments\nLine 1: Title. I’d suggest reviewing the title, as I think “main source of uncertainty” implies that you’ve compared biogeography to other sources of uncertainty in the land carbon cycle.\nLine 25: Worth making clear that you’ve reduced the spread of land carbon estimates for the present day, or does it also cover the future?\nLine 26: Again, emphasise that you’re talking about historical changes\nLine 33: the term “offsetting” tends to be used in carbon trading. I’d suggest being more quantitative here – perhaps cite the fraction of anthropogenic CO2 that the land absorbs?\nLine 54: Limitations could also include consistency of satellite timeseries, particularly as we’re interested in the response of the C cycle to increasing CO2 emissions\nLine 58-59: See line 1 comment.\nLine 76: Maybe clearer to say “global” GPP here\nLine 82: Would perhaps be better to assess how improved estimates of GPP affect the uncertainty in land carbon fluxes, and improve our ability to close the C budget over this period\nLine 87: Is the PFT total area also in km2 / year?\nLine 90-92: The premise of the whole paper hinges on this statement – see main reviewer comments.\nLine 93: Why didn’t the authors split EBF into tropical and temperature in a similar way to grasses? There would likely be a huge range in the GPP rate within the EBF PFT.\nLine 93-101: How do these values relate to the annual rate of GPP that would have come directly out of the models? Within a given PFT distribution, there would be a range of GPP rates (measured in similar units such as gC/year/sqm), so how are these rates different from just using the observed PFT distribution maps and taking a spatial mean and standard deviation of the modelled GPP rate?\nLine 101-103: This statement is probably for the discussion, and needs referencing / supporting with evidence from study\nLine 177: This relationship is bound to be strong if both variables co-vary strongly with area\nLine 229: Emphasise that it is PFT-specific total GPP\nLine 280-281: I don’t think you can conclude this from your results. I would argue that the higher uncertainty reduction in the approach presented in figs 1 & 2 is due to the over-dependence on area in the total GPP vs PFT area relationship\nLine 282-283: Again, I don’t think you can conclude this from your results. There is a lot of literature that shows parameter tuning from detailed site-specific observations - see for example Harper et al who used plant trait information to tune JULES, in places where the PFT fractions were very well documented. https://gmd.copernicus.org/articles/9/2415/2016/gmd-9-2415-2016.pdf\nLine 331: I think it’s important to be specific here - your study attempts to understand and reduce uncertainty in the *historical* to present-day carbon cycle. It highlights that DGVMs have work to do to improve the modelled spatial distribution of the their PFTs.\nLine 350-351: Worth acknowledging (apologies if I’ve missed it) that in using the S2 scenario, we would expect DGVMs to over-estimate the fraction of tropical broadleaf evergreen tree (which is the PFT class that normally has the highest contribution to global GPP)\nLine 359: See earlier comment on Tropical vs Temperate EBF\nLine 378-379: Maybe a better description of the metric is need here, given the >1200C criteria, I would assume it is rather something like accumulated temperature during the growing season?\nLine 417-418: Is this relationship on an annual timestep?\nLine 419: I think I’m right in saying that βi is the GPP rate per unit area here. In that case you would also need a sigma to sum across grid cells.", "reviewer_3": "Zhao et al present a method for constraining total global GPP and also to explore relationships between change in GPP over time with changes in PFT area. To constrain global GPP the authors define linear relationships across models between total global PFT specific GPP and PFT area and then constrain the range in global GPP values based on independent (though see caveat on that below) datasets of PFT area. In so doing, the authors highlight the importance of PFT area in modeled global GPP uncertainties. They also (for one independent dataset) constrain the mean global GPP to a higher value than is often reported, but which is more in line with other independent total GPP estimates.\nI served as a reviewer of the original version of this manuscript. I appreciate the extensive additional work the authors have done to revise the manuscript and to add extra, very involved analyses to address some of the concerns.\nIn my original review I emphasized that the goal of this study, to highlight the importance of differences in PFT area on uncertainties modeled GPP estimates, was much needed. I still think that, and furthermore I think the changes the authors have made to the paper have enhanced that message. I am very supportive of seeing this paper published.\nIn the original review I had concerns about the approach that have mostly been addressed. The difference in framing alone helps me to better digest and have greater confidence in the results and conclusions of this study. However, I still see some caveats that affect the strength of the conclusions drawn. I highlight these below, followed by key suggestions to address what I see as the remaining caveats. Following this general comments section I provide additional suggestions for several of the methods and results section before finishing with some minor comments.\nPlease forgive any repetition in the sections below. I have reviewed this paper over the course of the week because like this paper and I think it’s important and therefore I wanted to provide a thorough review. This required me thinking through some parts of the new analyses and responses to reviewers and took longer than I expected. I have tried to tidy up my comments but in a rush to get this review to you I may have missed some repetition. Apologies for that. I would be happy to discuss any of these comments further, either via the journal editors or directly via email.\n-- Natasha MacBean (nmacbean@uwo.ca)\nRemaining caveats with the approach used and therefore strength of the conclusions drawn:\n- I think the understanding gained from the linear relationships between total GPP and PFT area are more nuanced than the authors lay out in their paper (more details below).\n- I think the uncertainties on the PFT maps are too narrow and therefore the large reduction in total GPP spread obtained by using the PFT area from these two maps is not realistic (more details in methods and results below).\n- The fact that the two different LC maps result in widely different mean total GPP estimates, even with (what I think are) too narrow uncertainty bounds, highlights the issue that the reduction in spread of total GPP values from constraining with satellite derived PFT area (as shown in Figure 2) is unrealistically high. Therefore the strength put on this statement – the reduction of 86% – just does not seem to be reliable.\n- There is a circularity in using the ESA CCI map at least, because some models are using that map to prescribe PFTs.\nSuggestions to address these caveats:\n- I will let the authors digest what I have said about the R2 values (especially those not constrained to go through zero), the strength of the relationship between total GPP and PFT area, and what (I think) that says about our understanding of differences in PFT representation across models. But I think that there needs to be a more nuanced description and discussion of those results (and this section is currently very short). I would not place as much strength in the statement in lines 90-92 that “The strong relationships indicate that the total GPP of each PFT across DGVMs is largely dependent on its prescribed or simulated area, and that DGVMs simulate a similar per-area GPP for each PFT”.\n- I would like to see more realistic, wider uncertainties for each of the PFT products that included LC mapping uncertainties propagated through cross-walking uncertainties to the PFT maps (unless I’m missing this here – see more detailed comments in methods below), and therefore wider ranges of global total GPP (and a lower reduction on the total GPP spread) as a result. However, I would not wish that task of propagating cross-walking uncertainties to PFT maps on anyone! It was difficult enough in the Hartley et al. (2017) paper and we only explored some of the possible uncertainty “pathways” to demonstrate the issue – not the full range of uncertainties (which again leads me to think cross-walking uncertainties has not been accounted for here). So at the very least, I would like to see a clear statement and discussion that these error bars in Fig. 1 do not incorporate all sources of uncertainty in the final PFT maps, especially not in relation to the highly subjective cross-walking uncertainties, bioclimatic limits and somewhat arbitrary phenology thesholds used to link between LC class and PFT and to differentiate PFTs of different biomes, and therefore the real range could be much higher. Still, and this is a key point, even acknowledging wider uncertainties on the PFT maps, this could still constrain GPP to a higher value than is currently assumed (120-130PgC) – for MODIS at least.\n- Further to the previous point, I would like to see more discussion of all the sources of uncertainties in satellite derived PFT maps and what we as a community could do about them. This is mentioned in the discussion, but if there is space (I guess not) I would like to see more. Simply because the authors clearly do show that PFT area is an important source of uncertainty – that much is without a doubt for me. So which product is the best to use to prescribe PFTs in models, or to evaluate dynamic vegetation models (or to use as a constraint in the approach used in this study)? Do we just need an extensive worldwide campaign of validation of PFT maps – perhaps using citizen science? What else can be done? Examples of using using ancillary datasets of tree cover fraction, as has been done with ESA CCI, sound good to me, but I am sure there could be other approaches. Perhaps we could explore the idea of optical plant functional types (but see minor comment below in response to R1C3). Also we could use the approach you have taken with emulators to fully explore more extensive PFT map uncertainties.\nMethods – DGVMs\nIt’s not clear from the text, unless I missed it, that some of these models have prescribed, static PFT fractions and some have dynamic vegetation. It reads in the methods text as if all the models have dynamic vegetation, but that’s not the case is it? In fact, from reading the first section of the results it becomes clear that some models have prescribed and some simulate their dynamic vegetation. In this format of paper people may read the methods after the results, but for clarity I would add the information about whether the models have prescribed or simulated vegetation fractions to the DGVM methods section (and/or maybe Table S1)?\nFollowing on from this point – if most models have prescribed vegetation then the statement that parameters feedback on PFT distributions is not correct for those models that do not have dynamic vegetation enabled, and this confused the additional results that are presented in Figure S5. More on this below.\nMethods – emulators\nI really like the approach of developing PFT specific emulators to test PFT uncertainties. I appreciate this is an additional analysis prompted by a suggestion in my original review. I think it’s a very novel approach, so much so that there’s part of me that thinks it’s such a nice and novel analysis that it could be its own paper.\nHowever, I need more information here. What were the exact inputs to the emulator? Over how many grid cells was each PFT specific emulator developed (e.g. grid cells selected across a range of climate zones for each PFT?) and if so, how were those grid cells selected? Or was it a global emulator for each PFT? What GPP data were used to train the model? Over what time period was the outputs of the random forest model? Was a spin up and transient included in the emulator? Was each emulator independently tested/validated? If so, what were the results of that? All of this feeds into how we can interpret Fig. S5.\nMethods – Remote sensing estimates and uncertainties of global PFT distribution\nI appreciate the consideration the authors have taken in improving the methods section by adding details and clarifications. I am still a little uncertain how the uncertainties in both the underlying LC map, and the “cross-walking” to PFTs (including uncertainties in bioclimatic limits and phenology based thresholds to separate out different climate zones) has been considered. The authors discuss this, but a) I still don’t understand some details; and b) I still don’t think all these uncertainties have been considered.\nFor example:\n“We first obtained the uncertainties generated during mapping each land cover type by machine learning classifiers, and then assigned the uncertainties to PFTs based on translation rules from land cover to PFT maps (Bonan et al., 2002).”  \nWhat are the machine learning classifiers in this case? How does that give you the full uncertainty related to the LC map?\nAnd have uncertainties related to the translation rules between land cover and PFT map in the MODIS product also been considered – I can’t see any consideration of uncertainties in the cross-walking from LC to PFT in the Sulla-Menashe? In fact at Lines 407-408 the authors say “The uncertainty in MODIS PFTs was based on the reported uncertainties in the initial land cover mapping (Sulla-Menahshe et al.)”. This explicitly states that only uncertainties in the LC maps have been considered. Not any uncertainties in translating those (often broad) LC classes into PFT fractions. As we described in Hartley et al. (2017) the range of % you can give for some classes, especially in the very mixed and sparsely vegetated cover types, is almost impossible to know and very subjective. Similarly, you used additional bioclimatic limits from Bonan et al. (2002) and phenology based thresholds to further define PFTs according to their climate zone. However, these bioclimatic limits and phenology based (e.g. GDD) thresholds have uncertainties as well, and this is also not represented in your MODIS error bars. Some of the phenology based metrics – e.g., whether you cound GDD above a temperature threshold of 0 or 5°C – are very arbitrary.\nAlso:\n“Uncertainty for the ESA CCI product was derived by digitizing reported uncertainties”   what does digitizing reported uncertainties mean? So I see how the ancillary tree cover dataset can be used to estimate tree %, but what about all the other LC classes such as sparse vegetation – how did you assign fractions to each PFT? Did you use the ESA CCI cross-walking table for that? And therefore did you consider uncertainties in those very subjective cross-walking values, as described in Hartley et al., 2017.\nResults – Relationship between PFT GPP and area\nI understand perhaps better than I did in the original review what point the authors are trying to make here. Still, I would say that for some PFTs the linear relationship is not as strong (lower R2 values) and the fact that the R2 values are often lower when they are not forcing the regression line through zero also gives me pause. The logic for forcing the regression through zero is sound, but the authors are also just trying to see what the sensitivity of the model GPP is to the PFT area. In cases where the R2 is much lower, then actually what the non zero constrained regression is telling you is that actually the total GPP is not that sensitive to the area of the PFT as it is for other PFTs. This is especially the case for C3 grasses which are one of the most abundant PFTs worldwide (if you add up the area of the individual C3 categories in this study). So when the authors say in their response to R1C2 that “it is very surprising to see ALL models show EBF has similar productivity per area (i.e., the slopes of the total GPP – area relationships for each PFT in Fig. 1)” I would say I agree, but I also think that’s not the case for all PFTs.\nI think this analysis is interesting though because the differences across PFTs in terms of the strengths of the linear regression (unconstrained R2 values) provide insights into which PFTs have very similar representation across different models (and therefore where the prescribed PFT area is causing most of the difference in total area) versus the PFTs where area plays a role, but clearly other aspects of modeling that PFT are contributing to the spread in total GPP (which then gives lower R2 values).\nSome of the PFTs that clearly have a strong linear relationship between total GPP and area are not that surprising – i.e. C4 grasses have had little attention paid to them compared to C3 grasses, and so it is perhaps predictable that all models simulate a similar GPP per area as their underlying representation is very similar. In contrast, more work has been done on C3 photosynthesis – for example, different models for calculating C3 photosynthesis, different parameter values (some of which have been calibrated against more numerous datasets for C3 grasses), including optimality-based approaches, and including temperature acclimation, etc – and therefore I assume there are differences across models, and that explains the lower R2 values. (I also note that Figure S1 is not cited in the main text). Deciduous needleleaved forests are another example that have had relatively little attention paid to them compared to other PFTs and therefore the R2 values are quite high.\nResults – Updated land carbon cycle estimation using remotely sensed PFT area\nGiven the issues I mentioned above I am afraid I do not trust the narrow range of uncertainties used for either MODIS or ESA CCI. See my comments in the general section on suggestions for addressing the key caveats above.\nFig. 2: Caption include the fact that the DGVMs are from TRENDY and should say which version, even though this is in the methods (as many people won’t read the methods). In fact I would make it clear what is going into each of these estimates in the figure caption. Which products are going into MLEC for example, just FLUXCOM or a range of different upscaled EC products?\nResults – Contribution of PFT change to global GPP change\nThis section is really nice and interesting and yet it is a little lost in the overall message about constrained PFT area leading to constrained GPP based on the strength of the linear relationship between PFT GPP and area. I would have liked to see more here, although I appreciate there’s not much space in a letters journal. What is very curious about this is that we expect to see C4 grass outcompete at higher CO2 concentrations, if I’m remembering correctly. This is not directly discussed although the authors mention the potential water savings under elevated CO2. I would suggest the authors expand this though to suggest future studies check that that is indeed what is happening (water savings) and the incrased C3 is not because of incorrect process representation in the models.\nDiscussion\nLines 247-249: “To improve biogeography process in DGVMs, we suggest prioritizing efforts to clear hotspots of biogeography uncertainty (e.g., those identified in Fig. 3) and to develop a more reliable global PFT map”   Yes absolutely! Actually most of these models have no dynamic vegetation nor vegetation demography that is causing any kind of PFT change. They all have prescribed maps. So yes, a more reliable global PFT map is needed for models that do that! And yet, as the authors have shown, the two most widely available maps have very different areas, which, as the authors have also shown, feeds into differences in total GPP (though I just wouldn’t place as much emphasis on the narrowing of uncertainty in that total GPP estimate as the authors do). This is the key message of the paper here I think, but I believe it can be framed more clearly.\nLine 249-250: I guess I would say be careful with the use of “alternative” in this statement: “Using a consistent biogeography baseline, such as a PFT map from remote sensing as a diagnostic alternative” because these PFT maps are being used in models. I know ORCHIDEE and SDGVM use the same ESA CCI map (which also makes me wonder about the circularity in this approach – I would like the authors to comment on that).\nLines 256-261: This is only true for dynamic vegetation models, of which there are only 4 in the TRENDY v9 ensemble. Otherwise, photosynthesis rates etc do not feedback on PFT area. So I would revise this section accordingly.\nLines 269- : I like the additional emulator analysis and I really appreciate the authors attempting to address this comment from the first round of review. I have one minor comment and one concern:\nMy concern is that I don’t think the larger spread (smaller reduction in uncertainty in EM(DGVMS,MODIS) is because the models have tuned their structure to the wrong location. Or at least, I can’t see how that would be the case. Let’s say C4 grass is in the completely wrong location in many models because the bioclimatic limit used to separate out C3 and C4 grasses is wrong. If so, what does that imply, that their Vcmax parameter for example is adjusted to the climate of the wrong location to match in situ flux tower data? That seems convoluted to me, although I think maybe I see the argument. I’d be willing to bet though that most PFTs are not in the completely wrong location (mostly based on prescribed PFT maps). It would mean major biases in parameters or model structure to account for that. What I think is maybe more likely, or explains more of the difference between EM(DGVMs,MODIS) and DGVM_MOD is what I’ve said throughout this review – that the uncertainty reduction from constraining with MODIS (or ESA CCI) PFT map area is too high. I would guess that if you had more realistic uncertainties on your PFT maps, your constrined GPP would also be about 40-50% as in EM(DGVMs,MODIS). It still shows a big importance of PFT map, but again, I think it just cannot explain the majority of the spread across models. Still, this nice emulator result shows that differences simply due to thye PFT map are a greater issue than anyone would bet.\nMy minor comment is that I am not sure I would describe it as a difference between biogeography and biogeochemistry (Responses to R1C7 and R2C3 and new text describing Fig. S5) (with my issue being mostly with the use of biogeochemistry). If the location of the PFTs is wrong (e.g., Bias type 1 in Response to R2C3), then that’s also a biogeography issue right? Or am I misunderstanding? If the location is wrong then it could be that the bioclimatic envelopes for example are not right, and I would call that a biogeography “process” or limitation (that is intimately tied to area). On the other hand, you may have an issue with temperature optima, which would prohibit certain PFTs from growing in an area. I would call that a biogeochemistry (or plant physiology actually) process (could also be related to phenology thresholds or ecohydrology if not enough water was available). So it’s just a question of wording but I perhaps would more simply say that what you said in response to R2C3 – that there are two types of biases related to area and location, and that biases related to area only were accounted for by the analysis in Figure 1. But Figure 3 showed the second bias related to location distribution, and so you have done this analysis. Basically, I think the use of biogeochemistry here is confusing and also misleading given that for most models any “biogeochemistry/plant physiology” etc process is not feeding back on the distribution anyway. If the distirbution is wrong for most models then it’s because the prescribed PFT map is wrong (except for the dynamic vegetation models).\nReply to response to R1C3 in response to reviewers: I would perhaps disagree that PFTs can be directly linked to spectral signatures. In fact, it is very hard to distinguish trees from grasses, even in hyperspectral data and even when you add LiDAR data (as my PhD student Rubaya Pervin found in drylands given there’s also a strong spectral signature from dry grasses and/or biocrusts). In the past this has prompted the suggestion of using “optical functional types” (Ustin and Gamon, 2010) where global vegetation would be classified based on the clustering of spectral signatures. To my knowledge, no one has really tried this at regional or global scale. I would bet this would look different to the PFT maps we currently use and would be hard to match with anything we represent in models. So the question remains, how can we use satellite data to infer plant functional type fractional cover? The fusion of tree cover fraction in the ESA CCI product is a great advance. But it also depends on what is counted as a tree (height thresholds etc).\nReply to response to R1C6 in response to reviewers: Yes this is very interesting to see the range of estimates, and also the reporting of more recent studies that all have higher estimates. In fact, the higher mean value of GPP as constrained with MODIS is one of the most interesting things about this paper (although this is somewhat dampened by the fact that ESA CCI maintains the typical value of around 120PgC), although other methods have a higher value as the authors show in Fig. 2a. As I am signing this review I apologise for this shameless self-citation, but a paper in which I constrained GPP using SIF data (one of the earlier, simpler approach using just one model) also found a posterior (post parameter optimization) global value of just over 160PgC (MacBean et al., 2018). Please do not feel you need to include this is a reference, it is more for your interest. Since I wrote that paper in the back of my mind has been the thought that yes ok we constrained the high positive bias in ORCHIDEE (pre optimization it was above 190PgC), but we have not constrained it enough. I hadn’t quite processed until I read the authors response to the original review that there are indeed multiple studies pushing at this higher value (and hadn’t seen the Liu paper). As the authors say, the value of 120-130PgC has been talked about for Beer’s 2010 Science paper if not before and is the number that is most often quoted. It has been 130PgC in the the GCB for at least the last two years. So I think it’s great that the authors have moved away from a “correct” magnitude and instead show the uncertainties. Still, when you put all these studies together with the new study presented here, it does feel like the evidence is mounting for a higher value.\nMinor points:\nLine 29: would suggest “biogeography-climate” instead of “biogeography – climate” interactions. First time I read it I thought the words after the hyphen were making a different point.\nLines 54-56: another limitation for MODIS GPP (or any GPP derived from the production efficiency/LUE method) and the machine learning approaches I believe is the lack of accounting for the effect of photosynthesis on LUE (De Kauwe et al., 2016).\nLine 71: would maybe not use the word “prescribe” here (to be pedantic) – maybe just that in process-based models these processes are specific to each PFT ?\nI cannot see CLM5 in Fig. 1j? I’m guessing that’s why you have axes upper limits that are much higher than the main cluster of models\nI like the new additional analysis presented in Fig. S6. Nice!\nData Availability Statement needs to include where you got the TRENDY v9 GPP and PFT from as it’s actually not on the global carbon budget page that links to Mike O’Sullivans Github (https://mdosullivan.github.io/GCB/). Also, would you consider providing the emulators in a publicly available repository for people to test other PFT maps?\nReferences\nDe Kauwe, M. G., Keenan, T. F., Medlyn, B. E., Prentice, I. C. and Terrer. C. (2016) Satellite based estimates underestimate the effect of CO2 fertilisation on net primary productivity. Nature Climate Change, 6, 892-893.\nMacBean, N., Maignan, F., Bacour, C., Lewis, P., Peylin, P., Guanter, L., ... & Disney, M. (2018). Strong constraint on modelled global carbon uptake using solar-induced chlorophyll fluorescence data. Scientific reports, 8(1), 1973.\nUstin, S. L., & Gamon, J. A. (2010). Remote sensing of plant functional types. New Phytologist, 186(4), 795-816.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviewers\n============================================= =========================\nReviewer #1 (Remarks to the Author):\nR1C1: Land carbon sink remain largely uncertain among intercomparison models, thus restricting our capability of predicting, mitigating and adapting to climate change. There have been research efforts which leveraged the emergent constraint approach to reduce the predicted uncertainty for the variables (i.e., photosynthesis or gross/net primary production - GPP/NPP and ecosystem respiration – ER). The essence of emergent constraint is to identify the emergent relationships (typically linear) between a multi-model ensemble and an observational estimate. Here this study by Dr. Zhao et al have used this emergent constraint approach to identify the remarkably strong linear relationship between carbon flux (GPP and ER) and biogeographic area of vegetation biomes. Then the two sets of remote sensing datasets (MODIS and ESA CCI PFT maps) were used as the observational estimate to constraint the uncertainty of GPP and ER, respectively. I am highly impressive by the approach of using biogeographic vegetation biome area and its strong relationships because this offers significant insights in term of vegetation biome impact on carbon cycling.\nResponse: We thank the reviewer for their continued support of our study and for providing constructive feedback. We have endeavoured to use their comments to further improve our study.\nR1C2: Here I have some moderate suggestions to improve this study. First, the key uncertainty or weakness (if any) is that the remote sensing data of vegetation biome area has large uncertainty as shown in the results of figure 1 et al. The authors have made a good job to clarify this. Since the emergent constraint largely depends on the accuracy of observational variable – vegetation biome area here. The authors may consider clarifying that 1) the emergent constraint could still reduce the uncertainty even if the two sets of remote sensing datasets were used together by combing the remote sensing datasets (i.e., by getting the range of minimum and maximum values of GPP and ER); 2) the needs for future studies to improve or have more reliable observational estimate of vegetation biome area.\nResponse: Thank you for your suggestions. Please kindly see our detailed response below:\n1) We agree that the accuracy of remote sensing (RS) plant functional type (PFT) maps could affect the resulting GPP range. Though in the last round we have clarified that the uncertainty in RS PFT would not affect our main message as long as the strong GPP-area relationship stands, we now push further on this issue at the request of the reviewer. In this round, we collected five widely-used RS PFT maps (Table S3), including those from MODIS and ESA CCI.\nWe have now updated our main results (Fig. 1, 2, 3) using the ensemble of five RS PFT maps to constrain GPP instead of just using MODIS and ESA CCI PFT maps. After considering the inter-product difference in RS PFT maps, we found the GPP estimate constrained from 146.32 ± 26.39 PgC/year to 135.66 ± 6.54 PgC/year and ER from 143.39 ± 25.99 PgC/year to 133.01 ± 6.53 PgC/year, a 75% reduction for uncertainties in global GPP and ER estimates. This reduction is less than when only using the MODIS PFT map (86%), but is still very substantial and supports our main conclusion that vegetation biogeography is a main source of land carbon cycle modelling.\n2) We fully agree that it is critical for future studies to improve global PFT mapping, in particular future studies need to ensure the consistency across products to better support the land carbon cycle modelling. Using the newly collected five RS PFT maps, we identified the hotspots of biogeography disagreements among RS PFT maps (Fig. 3). We suggest that future study should prioritize efforts to clear these hotspots of PFT disagreements (L212-220) and our additional analysis also highlighted the PFTs that we need to work on, such as C3 cropland and shrubland (Table S3).\nTable S3. The PFT area estimates across five remote sensing-based products (unit: 106 km2). Only MODIS and ESA CCI PFT maps have reported part of their own uncertainty. σ refers to inter-product uncertainty (i.e., one standard deviation of the area from five products).\nESA CGLS_ GLC_ GLASS-\nPFTs MODIS unc unc unc unc unc σ\nCCI LC100 FCS30D GLC\nEBF 21.56 1.23 14.23 2.85 18.30 NA 15.30 NA 13.97 NA 3.23\nENFt 5.08 1.21 4.17 0.83 7.32 NA 5.24 NA 2.61 NA 1.72\nENFb 4.69 1.11 3.41 0.68 6.20 NA 5.13 NA 8.17 NA 1.79\nDBF 16.13 2.79 7.21 1.44 11.68 NA 16.57 NA 9.81 NA 4.04\nDNF 4.77 0.83 2.33 0.47 6.04 NA 6.32 NA 4.56 NA 1.58\nSHR 14.01 1.73 4.68 0.42 11.33 NA 14.17 NA 18.48 NA 5.08\nC3Gw 9.29 0.58 12.73 0.64 8.52 NA 9.96 NA 7.81 NA 1.90\nC3Gk 5.84 0.36 12.31 0.62 8.16 NA 6.46 NA 10.38 NA 2.71\nC3Gc 3.95 0.25 4.57 0.23 3.87 NA 5.08 NA 3.32 NA 0.68\nC3C 7.66 0.71 16.34 0.82 10.90 NA 13.61 NA 7.86 NA 3.74\nC4G 17.26 1.09 17.26 1.09 16.62 NA 17.26 NA 17.45 NA 0.32\nC4C 3.23 0.30 3.23 0.30 3.47 NA 3.23 NA 3.16 NA 0.12\nNote: the remote sensing-based products include MODIS v6.1 (Sulla-Menashe and Friedl, 2018), ESA CCI v2.0.8 (Harper et al., 2023), CGLS_LC100 v3 (Buchhorn et al., 2020), GLC_FCS30D (Zhang et al., 2024), and GLASS-GLC (Liu et al., 2020). Uncertainties (unc) for MODIS and ESA CCI PFT maps are derived from literature (see Methods). Uncertainties for other products are not available (NA).\nFigure 1. The relationships between total gross primary productivity (GPP) and total area for each PFT. The panels represent the relationships between GPP and area for (a) evergreen broadleaf forests (EBF), (b) temperate evergreen needleleaf forests (ENFt), (c) boreal evergreen needleleaf forests (ENFb), (d) deciduous broadleaf forests (DBF), (e) deciduous needleleaf forests (DNF), (f) shrubs (SHR), (g) warm C grasses (C3Gw), (h) cool 3 C grasses (C3Gk), (i) cold C grasses (C3Gc), (j) C crops (C3C), (k) C grasses (C4G), and 3 3 3 4 (l) C crops (C4C). The uncertainty of the relationship was quantified as the upper and lower 4 boundary of the relationship through a bootstrapping method for the regression (see Methods). The regression lines were forced through (0,0) as in theory there would be no GPP if PFT area was zero. The strength of the relationships was assessed by the coefficient of determination (R2). Meanwhile, we also provided R2 for the regression lines without forcing through (0,0) in parentheses. The vertical green lines stand for the average PFT areas from five remote sensing-based (RS) PFT maps, including MODIS v6.1 (Sulla-Menashe and Friedl, 2018), ESA CCI v2.0.8 (Harper et al., 2023), CGLS_LC100 v3 (Buchhorn et al., 2020), GLC_FCS30D (Zhang et al., 2024), and GLASS-GLC (Liu et al., 2020). The vertical shadings indicate the one standard deviation of the PFT area estimated by the five RS PFT maps. The horizontal green lines and shading intervals represent the constrained total GPP for each PFT using RS PFT areas.\nFigure 3. The difference in the distribution of PFTs from remote sensing-based (RS) products. The coefficient of variation (CV) for each PFT distribution across RS products is calculated as the ratio of the standard deviation to the mean PFT area, with higher CV values indicating greater uncertainties in PFT simulations. The abbreviations for PFTs and RS products correspond to those in Figure 1.\nR1C3: Second, the authors may also clarify the potential limitation of emergent constraint in term of underlying mechanisms underlying mechanisms (i.e., Winkler et al., 2019) and the spatial/temporal scales mismatch between data and models (i.e., Winkler et al., 2019). Ideally or necessarily, the research efforts in reducing the uncertainty of land carbon sink should be aimed by both improving the model structure or processes (i.e., the emergently vegetation biome distribution) and the emergent constraint approach here. This also related to the constraint of ER in this study because the current model structure have very simple structure for modeling ER which are strongly coupled with GPP.\nResponse: We hope to clarify that in the last round of revision, we have stepped away from stating our method as emergent constraint – we understand that there is a debate within the community now about the definition of emergent constraint (please see our replies to reviewers in the last round), and we did not intend to engage more into that discussion. Therefore, we now refer to our method as a simple correlation-based approach. The relationship we identified, strongly demonstrate the key role of vegetation biogeography on land carbon cycle modelling, and motivate future carbon cycle studies to consider biogeography (and their associated structures and parameters, L286-290) as a main place for improvements.\nR1C4: Third, the authors may also consider the potential effects of land usage impacts. Land use change likely have influenced the vegetation biome distribution more than climate change in term of spatial and temporal scales. Currently, it is not clear how the land usage change has been accounted for.\nResponse: Thank you. In this study, we focused on how inter-model difference in GPP modelling is influenced by inter-model difference in PFT distribution. We mainly used the S2 scenario estimates from TRENDY v9, which did not consider the effect of LUCC. Motivated by reviewer’s comment, we used the S3 scenario simulations, which include LUCC impact, to redo the analysis. We found the strong relationships between total GPP and total area for each PFT persist in S3 scenario (Fig. S9), which include the impact of LUCC. The new analysis indicates that the strong PFT GPP and Area relationships exist at the global scale in all scenarios. We have included a statement in L371-374 to clarify that.\nFigure S9. The relationships between total GPP and total area for each PFT. Here we presented the results using S3 scenario simulations from TRENDY v9.\nR1C5: Since the paper is very well written with insightful discussion, I also have the following minor suggestion for the authors to consider\nAbstract: may also consider reporting the results of ER\nResponse: Thank you for the positive comment on our discussion section. Following the suggestion, we have included the results about uncertainty reduced (i.e., approximately 75%) for ER using remote sensing-based PFT maps to further highlight the importance of biogeography for land carbon cycles (L25).\nR1C6: Lines 44-45: it would be important to mention competition, i.e., by accounting for size or cohort based competition for resources (light, nutrients and water)\nResponse: Thanks. We have added this point in L48.\n‘Considerable efforts have been made to identify the sources of uncertainty in models for land carbon cycle simulations, such as assessing climate forcing datasets and improving representations of ecological processes (e.g., thermal acclimation of photosynthesis, explicit representation of microbial processes, plant hydraulics, and competition for resources (Dantas De Paula et al., 2025)), and parameterizations (e.g., photosynthetic capacity and microbial carbon use efficiency) in process-based models.\nR1C7: Introduction: may consider clarifying what the emergent constraint is and the limitation of emergent constraint.\nResponse: Please see our response to R1C3. Since the last version, we have stepped away from referring to our method as emergent constraint as the request of other reviewers, to avoid potential confusions.\nR1C8: Figure 1: the color scheme would need to be improved. Also there are different definitions for vegetation biomes such as WWF vegetation biome. My knowledge is also that there is 1 km resolution of land cover maps https://www.earthenv.org/landcover I am not sure whether this is useful as an additional observational constraint.\nOverall, the whole manuscript has focused on the constraint of GPP, but the results of ER was indeed also important for carbon sink. Thus, the authors may consider balancing this.\nResponse: Thank you for the suggestions. 1) We have updated Fig. 1 (now Fig. S2) by adopting a more bright colour scheme, following the suggestion from the reviewer.\nFigure S2. The relationships between total GPP and total area for each PFT. Here we presented the results similar to Figure. 1, but showing the PFT areas observed from MODIS and ESA CCI only. The vertical shadings indicate the one standard deviation of each PFT map derived from existing literature (see Methods). The horizontal green and red lines and shading intervals represent the constrained total GPP for each PFT using MODIS or ESA CCI PFT area, and the uncertainty in GPP propagated from the PFT area.\n2) We agree with the reviewer on considering more RS PFT maps, such as the WWF vegetation biome classification and the land cover map that the reviewer mentioned. Please see our response to R1C2, in which we collected five widely-used global PFT maps, including the land cover that shared by the reviewer, to provide a more robust estimate of constrained GPP (Table S3 and Fig. 1).\nPlease kindly note that since WWF map does not adopt PFT classification and is not compatible with the PFT classifications in DGVMs, we ended up not using it. But we did include the one from the link that shared by the reviewer in our updated analysis.\n3) As for ER constraint, we agree with the reviewer that ER is important as GPP in determining the carbon sink. However, the challenge is that only a limited number of DGVMs provide PFT-specific ER; thus, we do not have sufficient data to do the analysis on ER at the same depth as GPP. As an alternative, we have constrained global ER based on the strong relationship between global GPP and global ER. We highlighted that this aspect needs improvement in the DGVM outputs (L394-397).\nR1C9: Figure 2: I am curious whether N deposition has been accounted for and how important it would be.\nResponse: That’s an interesting question – broadly the reviewer is asking about the impact of N deposition on vegetation biogeography change. Unfortunately, the current TRENDY DGVMs do not design enough scenarios to allow assessment of its impact on PFT distribution. However, considering that N deposition primarily occurs at a regional scale – particularly in regions with strong GHG and aerosol emissions, such as China and India (Zhu et al., 2025), it is likely to affect PFTs at a regional level. There is some evidence suggesting N deposition benefits certain species than others (Wang et al., 2021), however, we do not see any observational evidence at the biome level. We would like to regard this as an open question for regional level PFT change, and can potentially examine it once more data is available.\nR1C10: Discussion: as clarified, the authors may consider clarifying the limitation of emergent constraint and needs of improving model structure/processes for the mechanistically deriving/model the emergent properties of vegetation biome.\nThe section of discussion on structure (vegetation biome biogeographical distribution) vs biogeochemistry is interesting. It is also good to learn that the authors used emulator to disentangle their relative impacts. With this, the authors may consider elaborating a bit more for the implication of improving the model mechanisms in term of more research efforts on vegetation biome biogeographical distribution itself vs biogeochemistry.\nResponse: Please kindly refer to our response to R1C3 and R1C7, where we stated the reason for our withdrawal from using the term “emergent constraint” and clarified that we now refer to the relationship as a simple “GPP-area relationship”.\nWe appreciate the reviewer’s suggestion on elaborating more for the implication of improving model in term of biogeography vs biogeochemistry. Based on our emulator analysis, we have expanded our discussion regarding future efforts on these two aspects (L286-290).\n“Overall, the EM(DGVMs, MODIS) approach perhaps provides a more realistic and conservative estimate of biogeography's role in GPP modelling, and indicates that improvements in biogeography in GPP modelling should be accompanied by improvements in biogeochemistry.”\nReferences\nWinkler, A. J., Myneni, R. B. & Brovkin, V. Investigating the applicability of emergent constraints. Earth Syst. Dyn. 10, 501–523 (2019).\nReviewer #2 (Remarks to the Author):\nReviewer’s comments\nR2C1: This paper attempts to identify an observational constraint for dynamic global vegetation models (DGVMs) by establishing a linear regression relationship between plant function type (PFT) area and total gross primary productivity (GPP) for a given PFT. The authors use satellite-based observations to infer a plausible range of total global GPP per-PFT (GPPpft). They correctly identified that land surface models with dynamic vegetation components are often deficient when it comes to simulating present-day biogeography, and that using these PFT observations to constrain total global GPP estimates can tighten the range of GPP uncertainty.\nResponse: We appreciate the comments from the reviewer and have endeavoured to improve our manuscript following their suggestions. Please kindly see the point-by-point response below, where the line numbers refer to those in the clean version of the revised manuscript.\nR2C2: I have serious concerns, however that the linear regression on which this study is based is fundamentally flawed because both variables used in the regression have a strong dependence on PFT area. While I broadly agree with the conclusion that observed global area of a PFT can be used to refine global GPP estimates derived from DGVMs, I have serious reservations about the approach taken here. My criticism boils down to the fact that in this regression relationship, the independent variable (x-axis; PFT area) is also a central component of the dependent variable (y-axis; total global GPPpft), and therefore the resulting strong relationships shown for each PFT are not valid because it is essentially using PFT area to predict a variable that is strongly dependent on PFT area. In fact, the only reason why there is variation in the relationship (and not a perfect 1:1 line) is because there is variation in the per-model and per-grid cell GPPpft rate per unit area. In order to derive total global GPPpft, I assume the authors would have had to multiply the raw land surface model output GPPpft rate per unit area (likely in units similar to gC m-2 yr-1) by the grid cell PFT area, and then sum all grid cells (possibly also making a multi-year mean of the timeseries). Similarly, converting model output of PFT fraction to a total land area for a given PFT (Areapft) would involve multiplying the per grid cell fraction by each grid cell’s area (giving the exact value used in the GPP calculation in the previous sentence), and then summing all grid cells. If you take out the area dependence from the simple GPP calculation, you are left with the GPPpft rate again per unit area. Therefore, I think the linear regression on which the emergent constraint is based is flawed.\nResponse: We agree with the reviewer about the description of the calculation of PFT total area and total GPP, and we have used a similar approach in our analysis. We sensed, however, that there might be a misunderstanding regarding the interpretation of our results, which we clarify below.\nFor a given PFT in a given DGVM, the total GPP of the PFT can be simplified as:\nGPPpft = Aleaf*LAI*Area (1)\nwhere Aleaf is the average photosynthesis rate per leaf area (g/m2/yr), LAI is the leaf area index and Area is the total ground area for that PFT. Aleaf*LAI is the same as GPP rate that the reviewer mentioned.\nBased on equation (1), Area is certainly an important variable for GPP estimation (as noted by the reviewer). However, it does not mean a strong relationship between GPP and Area is a given. In fact, Aleaf and LAI variation between models are more often regarded as critical to GPP modelling in most DGVM uncertainty analysis, while the variation in Area is less noticed and considered to be not a main reason for inter-model variation in GPP. Our result directly presented that Area is a first-order factor controlling the GPPpft, and the impact of Area is stronger than Aleaf and LAI. We have provided the analysis in the last round of review to support the point (now Fig. S8a).\nFigure S8a. Comparison of the relationships between carbon cycle components and PFT characteristics. (a) The coefficient of determination (R2) for the relationships between PFT total GPP and PFT area (GPP & Area), or leaf area index (GPP & LAI) or leaf photosynthesis rate (GPP & Aleaf). The R2 is the coefficient of determination for the relationships, with higher values indicate a stronger correlation. The abbreviations for PFTs correspond to those in Figure 1.\nBy presenting the relationship between GPPpft and Area across models, our study offers a way to assess the importance of PFT area to the uncertainty (i.e., inter-model difference) in GPP modelling. We understand that the reviewer had another question on the variation in GPP rate within PFT which is more elaborated in R2C15. Please refer to our reply to R2C15 below.\nR2C3: Having said that, I do think there is value in highlighting that observations of PFT area can be used to constrain the simulated GPP of models in order to gain a more accurate measure of present-day GPP and consequently atmosphere-land carbon fluxes. I just don’t think that it can be done in the way that the authors have proposed. I believe that the authors would learn more about the models and the plausible constraints on GPP by looking in detail at the locations where PFT fractions match observations, and those where PFT fractions do not. I would also encourage the authors to examine other carbon-cycle emergent constraint studies and consider the definition presented in the Cox et al. (2019) paper:\n“The emergent constraint technique involves using the full ensemble of models to find an across-ensemble relationship between an observable feature of the Earth System (such as a trend, interannual variation or change in seasonality) and an uncertain aspect of the future.”\nI would argue that while PFT fraction is an observable feature of the Earth System (and GPP the uncertain aspect of the future), I don’t think that it provides a trend, interannual variation or change in seasonality, in the same way that other emergent constraint studies have. See for example, Cox et al. (2013) who found a strong linear relationship between the sensitivity of tropical land carbon storage to warming and the sensitivity of the annual growth rate of atmospheric CO2 to tropical temperature anomalies.\nResponse: Thank you for acknowledging the value in our study. We fully agree with the reviewer that there is a need to assess where the model PFT distribution is simulated differently. We have provided an analysis (updated Fig. S4, S5, S6) and a full discussion on that issue (L196-208). Also, we compared five widely-used PFT maps from remote sensing (RS), and quantified their disagreement over the world (updated Fig. 3). We have identified the key regions and key PFTs that needs improvement for mapping (L212-220).\nWe also hope to clarify that we have stepped away from referring to our method as “emergent constraint” in the last round, following the suggestion from the reviewer. We have now referred to our results as a “strong relationship between PFT GPP and PFT Area”.\nFigure 3. The difference in the distribution of PFTs from remote sensing-based (RS) products. The coefficient of variation (CV) for each PFT distribution across RS products is calculated as the ratio of the standard deviation to the mean PFT area, with higher CV values indicating greater uncertainties in PFT simulations. The abbreviations for PFTs and RS products correspond to those in Figure 1.\nOther Comments\nR2C4: Line 1: Title. I’d suggest reviewing the title, as I think “main source of uncertainty” implies that you’ve compared biogeography to other sources of uncertainty in the land carbon cycle.\nResponse: Thanks for your suggestion. Our full title conveys that “biogeography is a main source of uncertainty”. We understand biogeography is not “THE” main source of uncertainty since other biogeochemical factors are also important, but we believe our results are sufficient to indicate that the inter-model difference in PFT distribution substantially contributed to the inter-model difference in land carbon cycle modelling, and should be regarded as “a” main source of uncertainty.\nR2C5: Line 25: Worth making clear that you’ve reduced the spread of land carbon estimates for the present day, or does it also cover the future?\nResponse: Our analysis is based on the estimates of DGVMs in the past twenty years. To clarify, we have updated this sentence to ‘The dependence allowed us to reduce the spread of land carbon cycle estimates by approximately 75% using remote sensing-based PFT maps.’ We have stepped away from stating our approach as emergent constraint, given the term perhaps implies a future prediction.\nR2C6: Line 26: Again, emphasise that you’re talking about historical changes\nResponse: We have updated our statement to ‘We further found that 56 ± 21% of the climate-driven changes in global photosynthesis modelled by DGVMs were caused by changes in PFT distribution in the last two decades’, following the suggestion from the reviewer.\nR2C7: Line 33: the term “offsetting” tends to be used in carbon trading. I’d suggest being more quantitative here – perhaps cite the fraction of anthropogenic CO2 that the land absorbs?\nResponse: We have replaced ‘offsetting’ with ‘absorbing’ as suggested, and revised the sentence to ‘by absorbing approximately 25% of anthropogenic CO emissions (Friedlingstein 2 et al., 2023)’.\nR2C8: Line 54: Limitations could also include consistency of satellite timeseries, particularly as we’re interested in the response of the C cycle to increasing CO2 emissions\nResponse: We have added “lack of or inconsistent consideration of CO fertilization effect 2 (De Kauwe et al., 2016; Keenan et al., 2023)”, as suggested by the reviewer.\nR2C9: Line 58-59: See line 1 comment.\nResponse: Please kindly refer to our response to R2C2, where we clarified the key reasoning in our research.\nR2C10: Line 76: Maybe clearer to say “global” GPP here\nResponse: We have updated our statement to “We aim to quantify the dependence of PFT total GPP on PFT distribution and quantify the contribution of vegetation biogeography to the uncertainties in global GPP and ER estimation (see Methods)” for clarification.\nR2C11: Line 82: Would perhaps be better to assess how improved estimates of GPP affect the uncertainty in land carbon fluxes, and improve our ability to close the C budget over this period\nResponse: This is a really good point, and we did intend to conduct our analysis for the full C budget. However, a key limitation is that in TRENDY project output, many DGVMs only provide PFT-specific GPP but not ER or NEE, making it challenging to apply our approach to constrain ER or NEE for each PFT.\nAlternatively, we used constrained estimate of global GPP to narrow down global ER, based on the relationship between their global values across models (Fig. 2c). We ended up with a constrained global GPP of 135.66 ± 6.54 PgC/year and a constrained ER of 133.01 ± 6.53 PgC/year. We were unable to obtain a constrained NEE, as inter-model NEE is not correlated with inter-model GPP or ER. We have expanded our discussion (L394-397) to clarify this issue.\nR2C12: Line 87: Is the PFT total area also in km2 / year?\nResponse: The PFT total area reported here represents the average area over the period of 2001–2019, rather than an annual value or rate of change. Therefore, its unit is km², not km²/year. We hope to clarify that we are mostly interested in inter-model difference in this study. We have clarified this in L90.\nR2C13: Line 90-92: The premise of the whole paper hinges on this statement – see main reviewer comments.\nResponse: Please kindly refer to our response to R2C2 and R2C15.\nR2C14: Line 93: Why didn’t the authors split EBF into tropical and temperature in a similar way to grasses? There would likely be a huge range in the GPP rate within the EBF PFT.\nResponse: We did not split EBF into tropical and temperate EBF as most DGVMs regard EBF as one PFT and do not separate it into tropical and temperate, and most EBF is actually distributed in the tropics. We understand there would be a strong variation in GPP rate (g/m2/year) within one PFT, including EBF, however, analyzing such within-model variation is not the objective of our study. Our study focused on the inter-model discrepancy in GPP estimates, and therefore the inter-model difference in GPP rate and Area, not the spatial or temporal variations in GPP in one model.\nR2C15: Line 93-101: How do these values relate to the annual rate of GPP that would have come directly out of the models? Within a given PFT distribution, there would be a range of GPP rates (measured in similar units such as gC/year/sqm), so how are these rates different from just using the observed PFT distribution maps and taking a spatial mean and standard deviation of the modelled GPP rate?\nResponse: We apologize for the confusion – we noticed that a key point of misunderstanding is that we are discussing inter-model difference in GPP while the reviewer seemed to have concerns on the GPP rate variations within one PFT.\nUnfortunately, we cannot do “how are these rates different from just using the observed PFT distribution maps and taking a spatial mean and standard deviation of the modelled GPP rate” because DGVMs might not have simulations for a given PFT at places where remote sensing (RS) observations report there should exist that PFT. Even when there is overlap between DGVM PFT and RS PFT distributions, DGVM might simulate a very biased vegetation status (such as LAI) for a given PFT in that region, and consequently result in a biased GPP rate. Therefore, it would be very challenging to conduct the analysis from the pixel level (pixel level GPP rate for each PFT). Instead, our analysis is conducted at the biome scale (total GPP and total area for each PFT).\nTo avoid misunderstanding of the results, we have rephrased the statement to ‘For each of the twelve PFT types on the land surface, we observed a strong and significant relationship (p < 0.001) between GPP and area across models (Fig. 1). The strong relationships indicate that the inter-DGVMs discrepancy of total GPP of each PFT is largely caused by the difference in the prescribed or simulated PFT area of DGVMs, and that DGVMs simulate a similar per-area GPP on average for each PFT (i.e., the slopes of the GPP-area relationships in Fig. 1).’\nR2C16: Line 101-103: This statement is probably for the discussion, and needs referencing / supporting with evidence from study\nResponse: We have now moved the statement to discussion (L307-309) and added references (Hu et al., 2018; Rogers et al., 2017), following the suggestion from the reviewer.\nR2C17: Line 177: This relationship is bound to be strong if both variables co-vary strongly with area\nResponse: We appreciate your comment; however, we believe that the strong correlation between global GPP and global ER is not merely driven by their covariance with area. As we clarified in our response to R2C2, the strong relationship between PFT total GPP and total area across DGVMs is not a given, nor is the relationship between ER and area. The observed strong relationship between global GPP and global ER, indicates a strong ecological linkage between carbon uptake (GPP) and carbon release (ER), though in reality the relationship might not be that strong, considering the different factors controlling GPP and ER. The strong relationship between global GPP and global ER across models is likely due to a carbon sink centric view adopted by DGVMs, where the input of the carbon through photosynthesis largely determined the magnitude of ER (Fatichi et al., 2019), particularly over the long term.\nR2C18: Line 229: Emphasise that it is PFT-specific total GPP\nResponse: We have updated the statement to “PFT-specific total GPP” in L232.\nR2C19: Line 280-281: I don’t think you can conclude this from your results. I would argue that the higher uncertainty reduction in the approach presented in figs 1 & 2 is due to the over-dependence on area in the total GPP vs PFT area relationship\nResponse: Please kindly refer our response to R2C2 regarding the dependence of PFT total GPP on PFT total area.\nR2C20: Line 282-283: Again, I don’t think you can conclude this from your results. There is a lot of literature that shows parameter tuning from detailed site-specific observations - see for example Harper et al who used plant trait information to tune JULES, in places where the PFT fractions were very well documented. https://gmd.copernicus.org/articles/9/2415/2016/gmd-9-2415-2016.pdf\nResponse: Thank you for sharing the study with us. We have removed the statement as we found it challenging to make a generalised statement on the strategies of tuning parameters for different DGVMs.\nR2C21: Line 331: I think it’s important to be specific here - your study attempts to understand and reduce uncertainty in the *historical* to present-day carbon cycle. It highlights that DGVMs have work to do to improve the modelled spatial distribution of the their PFTs.\nResponse: Thank you for your suggestion. We have stated that we found the inter-model difference in land carbon cycle simulation is related to inter-model difference in PFT distribution, and that issue should persist regardless of whether we use DGVM simulations for past, present or future, though our analysis used DGVM simulations for the past twenty years. We have stated in our method we used TRENDY v9 simulations from 2000 to 2019 to be compatible with remote sensing-based PFT maps.\nR2C22: Line 350-351: Worth acknowledging (apologies if I’ve missed it) that in using the S2 scenario, we would expect DGVMs to over-estimate the fraction of tropical broadleaf evergreen tree (which is the PFT class that normally has the highest contribution to global GPP)\nResponse: Thanks for pointing this out. We have revised this sentence to acknowledge we used S2 scenario in our main analysis. Moreover, by doing additional analysis using S3 scenario, which further includes the impact of land use and land cover change, we have confirmed that our main findings remain robust (Fig. S9). We clarified this point in L370-374, “We mainly used the simulations from the S2 scenario to examine the relationship between PFT-specific total GPP and area. Changing the simulations to other scenarios, such as S3, would not affect the relationship between PFT total GPP and PFT area across models (Fig. S9), however we would like to exclude the consideration of land use and land cover change in this analysis.”\nFigure S9. The relationships between total GPP and total area for each PFT. Here we presented the results using S3 scenario simulations from TRENDY v9.\nR2C23: Line 359: See earlier comment on Tropical vs Temperate EBF\nResponse: Please kindly refer to our response to R2C14.\nR2C24: Line 378-379: Maybe a better description of the metric is need here, given the >1200C criteria, I would assume it is rather something like accumulated temperature during the growing season?\nResponse: We apologize for the unclear description here. Based on Prentice et al (1992), we have revised the statement to ‘an accumulated growing degree-days above 5°C (GDD ) 5 exceeding 1200°C’.\nR2C25: Line 417-418: Is this relationship on an annual timestep?\nResponse: We examined this relationship based on multi-year (2001-2019) average GPP pft and Area . We have clarified this by revising this sentence to ‘We obtained the multi-year pft average total GPP and total area for each PFT over the period 2001-2019, and examined the relationships between them’.\nR2C26: Line 419: I think I’m right in saying that βi is the GPP rate per unit area here. In that case you would also need a sigma to sum across grid cells.\nResponse: As clarified in our response to R2C2, βi represents the dependence of total GPP on total area for ith PFT (i = 1, 2, 3, …, 12, representing EBF, ENFt, ENFb, ..., C3C) across models (i.e., the slopes of the total GPP – total area relationships for each PFT in Fig. 1). We already summed pixel level values to the biome scale when calculating PFT total GPP (unit: PgC/year) and PFT total area (km2), separately. Therefore, an additional summation (sigma) across grid cells is not required. Please also kindly refer to our answers to R2C2 and R2C15, where we hope to clarify this point of confusion in the review process.\nReviewer #3 (Remarks to the Author):\nR3C1: Zhao et al present a method for constraining total global GPP and also to explore relationships between change in GPP over time with changes in PFT area. To constrain global GPP the authors define linear relationships across models between total global PFT specific GPP and PFT area and then constrain the range in global GPP values based on independent (though see caveat on that below) datasets of PFT area. In so doing, the authors highlight the importance of PFT area in modeled global GPP uncertainties. They also (for one independent dataset) constrain the mean global GPP to a higher value than is often reported, but which is more in line with other independent total GPP estimates.\nI served as a reviewer of the original version of this manuscript. I appreciate the extensive additional work the authors have done to revise the manuscript and to add extra, very involved analyses to address some of the concerns.\nIn my original review I emphasized that the goal of this study, to highlight the importance of differences in PFT area on uncertainties modeled GPP estimates, was much needed. I still think that, and furthermore I think the changes the authors have made to the paper have enhanced that message. I am very supportive of seeing this paper published.\nIn the original review I had concerns about the approach that have mostly been addressed. The difference in framing alone helps me to better digest and have greater confidence in the results and conclusions of this study. However, I still see some caveats that affect the strength of the conclusions drawn. I highlight these below, followed by key suggestions to address what I see as the remaining caveats. Following this general comments section I provide additional suggestions for several of the methods and results section before finishing with some minor comments.\nPlease forgive any repetition in the sections below. I have reviewed this paper over the course of the week because like this paper and I think it’s important and therefore I wanted to provide a thorough review. This required me thinking through some parts of the new analyses and responses to reviewers and took longer than I expected. I have tried to tidy up my comments but in a rush to get this review to you I may have missed some repetition. Apologies for that. I would be happy to discuss any of these comments further, either via the journal editors or directly via email.\n-- Natasha MacBean (nmacbean@uwo.ca)\nResponse: We thank the reviewers 3 and 4 very much for their time and effort in reviewing our manuscript and for their thoughtful comments. We are glad to learn that we have resolved most of their concerns and they are supportive of seeing this paper published. The reviewers provide some very in-depth thoughts, and we would be happy to clarify, and in some cases, explore with them below.\nR3C2: Remaining caveats with the approach used and therefore strength of the conclusions drawn:\n- R3C2.1: I think the understanding gained from the linear relationships between total GPP and PFT area are more nuanced than the authors lay out in their paper (more details below).\n- R3C2.2: I think the uncertainties on the PFT maps are too narrow and therefore the large reduction in total GPP spread obtained by using the PFT area from these two maps is not realistic (more details in methods and results below).\n- R3C2.3: The fact that the two different LC maps result in widely different mean total GPP estimates, even with (what I think are) too narrow uncertainty bounds, highlights the issue that the reduction in spread of total GPP values from constraining with satellite derived PFT area (as shown in Figure 2) is unrealistically high. Therefore the strength put on this statement – the reduction of 86% – just does not seem to be reliable.\n- R3C2.4: There is a circularity in using the ESA CCI map at least, because some models are using that map to prescribe PFTs.\nSuggestions to address these caveats:\n- I will let the authors digest what I have said about the R2 values (especially those not constrained to go through zero), the strength of the relationship between total GPP and PFT area, and what (I think) that says about our understanding of differences in PFT representation across models. But I think that there needs to be a more nuanced description and discussion of those results (and this section is currently very short). I would not place as much strength in the statement in lines 90-92 that “The strong relationships indicate that the total GPP of each PFT across DGVMs is largely dependent on its prescribed or simulated area, and that DGVMs simulate a similar per-area GPP for each PFT”.\n- I would like to see more realistic, wider uncertainties for each of the PFT products that included LC mapping uncertainties propagated through cross-walking uncertainties to the PFT maps (unless I’m missing this here – see more detailed comments in methods below), and therefore wider ranges of global total GPP (and a lower reduction on the total GPP spread) as a result. However, I would not wish that task of propagating cross-walking uncertainties to PFT maps on anyone! It was difficult enough in the Hartley et al. (2017) paper and we only explored some of the possible uncertainty “pathways” to demonstrate the issue – not the full range of uncertainties (which again leads me to think cross-walking uncertainties has not been accounted for here). So at the very least, I would like to see a clear statement and discussion that these error bars in Fig. 1 do not incorporate all sources of uncertainty in the final PFT maps, especially not in relation to the highly subjective cross-walking uncertainties, bioclimatic limits and somewhat arbitrary phenology thesholds used to link between LC class and PFT and to differentiate PFTs of different biomes, and therefore the real range could be much higher. Still, and this is a key point, even acknowledging wider uncertainties on the PFT maps, this could still constrain GPP to a higher value than is currently assumed (120-130PgC) – for MODIS at least.\n- Further to the previous point, I would like to see more discussion of all the sources of uncertainties in satellite derived PFT maps and what we as a community could do about them. This is mentioned in the discussion, but if there is space (I guess not) I would like to see more. Simply because the authors clearly do show that PFT area is an important source of uncertainty – that much is without a doubt for me. So which product is the best to use to prescribe PFTs in models, or to evaluate dynamic vegetation models (or to use as a constraint in the approach used in this study)? Do we just need an extensive worldwide campaign of validation of PFT maps – perhaps using citizen science? What else can be done? Examples of using using ancillary datasets of tree cover fraction, as has been done with ESA CCI, sound good to me, but I am sure there could be other approaches. Perhaps we could explore the idea of optical plant functional types (but see minor comment below in response to R1C3). Also we could use the approach you have taken with emulators to fully explore more extensive PFT map uncertainties.\nResponse: We appreciate the comments from the reviewer and their suggestions on how to address the comments. Following their suggestions, we have done the following analyses.\nR3C2.1 – We agree that the relationship of total PFT area and total GPP is not always strong. We therefore expanded our discussion in L303-307 to clarify that, and highlight some PFTs that have relatively weak correlations (such as C grasslands). For those PFTs without a 3 strong relationship, we expect their GPP uncertainty is more contributed by the biogeochemical process in DGVMs.\nR3C2.2 – We fully agree with the reviewer regarding the uncertainty in remote sensing-based (RS) PFT maps, it is somewhat disheartening to see that after a few decades of development, we still cannot get a reliable/consistent RS PFT map. The uncertainty in RS PFT maps comes either from their own production process (e.g., cross-walking) or from the difference across products in their classification criteria. To fully assess the uncertainty, we now collected five widely-used RS PFT maps and listed their own and inter-product uncertainty. We are now using inter-product PFT uncertainty in the main analysis since many products actually do not have “own production uncertainty” reported (Table S3). We added a full discussion on the issue (L334-344).\nAfter some consideration, we have decided to update our main results by using the mean of five RS PFT maps and the inter-product uncertainty in Figure 1. Using the ensemble of RS PFT maps, we still find a 75% decrease in GPP uncertainty. Though this reduction is less than when only using the MODIS PFT and ESA CCI maps (86%), but is still very substantial and supports our main conclusion that vegetation biogeography is a main source of land carbon cycle modelling.\nWe appreciate the reviewer’s understanding that it is beyond our effort to discuss RS PFTs mapping in this study, as our key is to demonstrate the strong role of PFT area in influencing carbon cycle modelling, which is largely neglected by previous studies. We do intend to further explore PFT mapping in future and would like to discuss with the reviewers through other channels.\nTable S3. The PFT area estimates across five remote sensing-based products (unit: 106 km2). Only MODIS and ESA CCI PFT maps have reported part of their own uncertainty. σ refers to inter-product uncertainty (i.e., one standard deviation of the area from five products).\nESA CGLS_ GLC_ GLASS-\nPFTs MODIS unc unc unc unc unc σ\nCCI LC100 FCS30D GLC\nEBF 21.56 1.23 14.23 2.85 18.30 NA 15.30 NA 13.97 NA 3.23\nENFt 5.08 1.21 4.17 0.83 7.32 NA 5.24 NA 2.61 NA 1.72\nENFb 4.69 1.11 3.41 0.68 6.20 NA 5.13 NA 8.17 NA 1.79\nDBF 16.13 2.79 7.21 1.44 11.68 NA 16.57 NA 9.81 NA 4.04\nDNF 4.77 0.83 2.33 0.47 6.04 NA 6.32 NA 4.56 NA 1.58\nSHR 14.01 1.73 4.68 0.42 11.33 NA 14.17 NA 18.48 NA 5.08\nC3Gw 9.29 0.58 12.73 0.64 8.52 NA 9.96 NA 7.81 NA 1.90\nC3Gk 5.84 0.36 12.31 0.62 8.16 NA 6.46 NA 10.38 NA 2.71\nC3Gc 3.95 0.25 4.57 0.23 3.87 NA 5.08 NA 3.32 NA 0.68\nC3C 7.66 0.71 16.34 0.82 10.90 NA 13.61 NA 7.86 NA 3.74\nC4G 17.26 1.09 17.26 1.09 16.62 NA 17.26 NA 17.45 NA 0.32\nC4C 3.23 0.30 3.23 0.30 3.47 NA 3.23 NA 3.16 NA 0.12\nNote: the remote sensing-based products include MODIS v6.1 (Sulla-Menashe and Friedl, 2018), ESA CCI v2.0.8 (Harper et al., 2023), CGLS_LC100 v3 (Buchhorn et al., 2020), GLC_FCS30D (Zhang et al., 2024), and GLASS-GLC (Liu et al., 2020). Uncertainties (unc) for MODIS and ESA CCI PFT maps are derived from literature (see Methods). Uncertainties for other products are not available (NA).\nFigure 1. The relationships between total gross primary productivity (GPP) and total area for each PFT. The panels represent the relationships between GPP and area for (a) evergreen broadleaf forests (EBF), (b) temperate evergreen needleleaf forests (ENFt), (c) boreal evergreen needleleaf forests (ENFb), (d) deciduous broadleaf forests (DBF), (e) deciduous needleleaf forests (DNF), (f) shrubs (SHR), (g) warm C grasses (C3Gw), (h) cool 3 C grasses (C3Gk), (i) cold C grasses (C3Gc), (j) C crops (C3C), (k) C grasses (C4G), and 3 3 3 4 (l) C crops (C4C). The uncertainty of the relationship was quantified as the upper and lower 4 boundary of the relationship through a bootstrapping method for the regression (see Methods). The regression lines were forced through (0,0) as in theory there would be no GPP if PFT area was zero. The strength of the relationships was assessed by the coefficient of determination (R2). Meanwhile, we also provided R2 for the regression lines without forcing through (0,0) in parentheses. The vertical green lines stand for the average PFT areas from five remote sensing-based (RS) PFT maps, including MODIS v6.1 (Sulla-Menashe and Friedl, 2018), ESA CCI v2.0.8 (Harper et al., 2023), CGLS_LC100 v3 (Buchhorn et al., 2020), GLC_FCS30D (Zhang et al., 2024), and GLASS-GLC (Liu et al., 2020). The vertical shadings indicate the one standard deviation of the PFT area estimated by the five RS PFT maps. The horizontal green lines and shading intervals represent the constrained total GPP for each PFT using RS PFT areas.\nR3C2.3 – Please kindly refer to our response to R3C2.2. We have now included inter-product uncertainty for RS PFT maps and reported a more conservative estimate of reduction in uncertainty. The GPP uncertainty has now reduced from 26.39 PgC/year to 6.54 PgC/year.\nR3C2.4 – Please kindly refer to our response to R3C2.2. Though appraising RS PFT maps is beyond our intention, our results indicate the place where we can further improve PFT mapping. From Fig. 1, we observed a large spread of C3C. From the updated Fig. 3, we identified the hotspots for PFT disagreements between the five RS PFT maps, including EBF in eastern U.S. and western Europe, ENF in northern Canada and eastern Russia, DNF in North America and Europe, and SHR and C3G across the globe. It would be beneficial to collect more samples there for future mapping and investigation.\nFigure 3. The difference in the distribution of PFTs from remote sensing-based (RS) products. The coefficient of variation (CV) for each PFT distribution across RS products is calculated as the ratio of the standard deviation to the mean PFT area, with higher CV values indicating greater uncertainties in PFT simulations. The abbreviations for PFTs and RS products correspond to those in Figure 1.\nR3C3: Methods – DGVMs\nIt’s not clear from the text, unless I missed it, that some of these models have prescribed, static PFT fractions and some have dynamic vegetation. It reads in the methods text as if all the models have dynamic vegetation, but that’s not the case is it? In fact, from reading the first section of the results it becomes clear that some models have prescribed and some simulate their dynamic vegetation. In this format of paper people may read the methods after the results, but for clarity I would add the information about whether the models have prescribed or simulated vegetation fractions to the DGVM methods section (and/or maybe Table S1)?\nFollowing on from this point – if most models have prescribed vegetation then the statement that parameters feedback on PFT distributions is not correct for those models that do not have dynamic vegetation enabled, and this confused the additional results that are presented in Figure S5. More on this below.\nResponse: We have clarified in Method section and updated Table S1 to indicate whether DGVMs enabled dynamic biogeography modelling or not. In short, most DGVMs prescribed static PFT fractions, except for JULES, LPJ-GUESS, LPX-Bern, and SDGVM.\nTable S1. The dynamic global vegetation models (DGVMs) from TRENDY v9, and the availability of their simulations of the area and gross primary productivity (GPP) for different plant functional types (PFTs) a.\nPFTs\nPFT\nDGVMs Note\ndistribution\nEBF ENF DBF DNF SHR C3G C4G C3C C4C\nCABLE YES YES YES YES YES YES YES NO NO Prescribed\nCLASS YES YES YES YES NO YES YES YES YES Prescribed\nC4C includes Corn, Millet,\nCLM5 YES YES YES YES YES YES YES YES YES Prescribed\nSorghum and Sugarcane\nC3C includes C3 crop and\nJSBACH YES YES YES YES YES YES YES YES YES Prescribed C3 pasture; C4C includes\nC4 crop and pasture\nC3C and C4C area the\nJULES YES YES YES YES YES YES YES YES YES Simulated\nsame as JSBACH\nDBF includes boreal\nneedleaf summergreen,\nLPJ_ temperate broadleaf\nYES YES YES YES NO YES YES YES YES Simulated\nGUESS summergreen, temperate\nshade-intolerant broadleaf\nsummergreen\nDBF includes tropical\nbroad raingreen,\ntemperate and boreal\nLPX_Bern YES YES YES YES NO YES YES YES YES Simulated broad summergreen C3C includes cropland C3\nherbaceous and pasture\nC3 herbaceous; C4C\nincludes cropland C4\nherbaceous and pasture\nC4 herbaceous\nDBF includes tropical\nbroadleaf raingreen,\nORCHIDEE YES YES YES YES NO YES YES YES YES Prescribed\ntemperate and boreal\nbroadleaf summergreen\nSame as ORCHIDEE, but\nORCHIDEE\nYES YES YES YES NO YES YES YES YES Prescribed not used, as GPP values\n-CNP\nare quite low\nSDGVM YES YES YES YES NO YES YES YES YES Simulated\nYIBs YES YES YES NO YES YES YES YES YES Prescribed\na The PFTs include evergreen needleleaf forests (ENF), evergreen broadleaf forests (EBF), deciduous needleleaf forests (DNF), deciduous broadleaf forests (DBF), shrubs (SHR), C 3 grasses (C3G), C grasses (C4G), C crops (C3C), and C crops (C4C).\n4 3 4\nIt is a great point from the reviewer that, for models that do not have dynamic vegetation enabled, the parameters feedback on PFT distribution is not possible. We have updated the statement accordingly to reflect this point (L258). However, since our main result is about inter-model GPP difference caused by inter-model difference in PFT distribution, we are less focused on exploring the dynamic changes in PFT over time for DGVMs. Therefore, in figure S5, we are simply discussing a case if we replace PFTs (even static) in DGVMs with a consistent RS PFT, what would the estimated GPP be.\nR3C4: Methods – emulators\nI really like the approach of developing PFT specific emulators to test PFT uncertainties. I appreciate this is an additional analysis prompted by a suggestion in my original review. I think it’s a very novel approach, so much so that there’s part of me that thinks it’s such a nice and novel analysis that it could be its own paper.\nHowever, I need more information here. What were the exact inputs to the emulator? Over how many grid cells was each PFT specific emulator developed (e.g. grid cells selected across a range of climate zones for each PFT?) and if so, how were those grid cells selected? Or was it a global emulator for each PFT? What GPP data were used to train the model? Over what time period was the outputs of the random forest model? Was a spin up and transient included in the emulator? Was each emulator independently tested/validated? If so, what were the results of that? All of this feeds into how we can interpret Fig. S5.\nResponse: We also feel excited about emulators, as they seem to suggest we can at least use machine learning to quickly “run” DGVMs. We are happy to provide more details below:\nSpecifically, the inputs to each emulator included annual climate forcings (temperature, precipitation, radiation, pressure and wind speed) during 2001-2019 from CRU JRA v2.1 datasets, which was consistent with the climate forcings for DGVMs in TRENDY v9. Each emulator was constructed globally for each PFT, using pixel PFT GPP outputs derived from the respective DGVMs in each year. We excluded the LPJ-GUESS model from our emulator analysis, as it reports GPP rates in per unit of grid cell, rather than per unit of PFT area like other models. We did not explicitly include a spin-up period, as the emulator training was directly based on steady-state annual DGVM outputs. Additionally, each emulator was independently validated through cross-validation, showing consistently high predictive accuracy (mean R² values > 0.8 for all PFTs, Table S4). We have now included Table S4 and methodological details in the revised manuscript (L480-493).\nTable S4. The coefficient of determination (R2) of emulators for PFT GPP simulations\nLPX- ORCHID\nCABLE CLASSIC CLM5 JSBACH JULES SDGVM YIBs\nBern EE\nEBF cal 0.98 0.99 0.87 0.99 0.98 0.38 0.93 0.99 0.95\nval 0.98 0.99 0.86 0.99 0.98 0.54 0.93 0.98 0.95\nENFt cal 0.97 0.98 0.99 0.89 0.94 0.27 0.85 0.97 0.74\nval 0.97 0.98 0.99 0.88 0.94 0.76 0.85 0.97 0.72\nENFb cal 0.96 0.98 0.98 0.80 0.93 0.95 0.84 0.89 0.83\nval 0.96 0.97 0.97 0.78 0.91 0.95 0.84 0.88 0.82\nDBF cal 0.97 0.99 0.97 0.99 0.98 0.46 0.91 0.95 0.77\nval 0.97 0.99 0.97 0.99 0.98 0.84 0.91 0.95 0.76\nDNF cal 0.98 0.97 0.92 0.87 0.94 0.91 0.61 0.93 0.77\nval 0.98 0.97 0.92 0.85 0.94 0.86 0.61 0.93 0.76\nSHR cal 0.98 0.98 0.95 0.99 0.97 0.44 0.84 0.96 0.66\nval 0.98 0.98 0.95 0.99 0.97 0.38 0.84 0.96 0.66\nC3Gw cal 0.96 1.00 0.96 0.98 0.99 0.77 0.78 0.90 0.70\nval 0.96 0.99 0.96 0.97 0.99 0.60 0.78 0.90 0.69\nC3Gk cal 0.96 0.98 0.91 0.94 0.98 0.44 0.84 0.96 0.59\nval 0.96 0.98 0.91 0.94 0.98 0.38 0.84 0.96 0.59\nC3Gc cal 0.97 0.99 0.91 0.96 0.96 0.86 0.65 0.92 0.80\nval 0.94 0.98 0.91 0.95 0.95 0.86 0.64 0.91 0.77\nC3C cal 0.96 0.96 0.98 0.97 0.99 0.88 0.84 0.95 0.74\nval 0.96 0.95 0.97 0.97 0.99 0.84 0.83 0.95 0.73\nC4G cal 0.97 0.99 0.91 0.96 1.00 0.93 0.91 0.84 0.77\nval 0.97 0.99 0.90 0.96 1.00 0.92 0.91 0.83 0.77\nC4C cal 0.97 0.95 0.90 0.95 0.99 0.84 0.86 0.96 0.80\nval 0.97 0.94 0.90 0.94 0.99 0.94 0.86 0.96 0.78\nNote: cal and val indicate the R2 for training and validation, respectively.\nR3C5: Methods – Remote sensing estimates and uncertainties of global PFT distribution\nI appreciate the consideration the authors have taken in improving the methods section by adding details and clarifications. I am still a little uncertain how the uncertainties in both the underlying LC map, and the “cross-walking” to PFTs (including uncertainties in bioclimatic limits and phenology based thresholds to separate out different climate zones) has been considered. The authors discuss this, but a) I still don’t understand some details; and b) I still don’t think all these uncertainties have been considered.\nFor example:\n“We first obtained the uncertainties generated during mapping each land cover type by machine learning classifiers, and then assigned the uncertainties to PFTs based on translation rules from land cover to PFT maps (Bonan et al., 2002).” à What are the machine learning classifiers in this case? How does that give you the full uncertainty related to the LC map? And have uncertainties related to the translation rules between land cover and PFT map in the MODIS product also been considered – I can’t see any consideration of uncertainties in the cross-walking from LC to PFT in the Sulla-Menashe? In fact at Lines 407-408 the authors say “The uncertainty in MODIS PFTs was based on the reported uncertainties in the initial land cover mapping (Sulla-Menahshe et al.)”. This explicitly states that only uncertainties in the LC maps have been considered. Not any uncertainties in translating those (often broad) LC classes into PFT fractions. As we described in Hartley et al. (2017) the range of % you can give for some classes, especially in the very mixed and sparsely vegetated cover types, is almost impossible to know and very subjective. Similarly, you used additional bioclimatic limits from Bonan et al. (2002) and phenology based thresholds to further define PFTs according to their climate zone. However, these bioclimatic limits and phenology based (e.g. GDD) thresholds have uncertainties as well, and this is also not represented in your MODIS error bars. Some of the phenology based metrics – e.g., whether you cound GDD above a temperature threshold of 0 or 5°C – are very arbitrary.\nAlso:\n“Uncertainty for the ESA CCI product was derived by digitizing reported uncertainties” à what does digitizing reported uncertainties mean? So I see how the ancillary tree cover dataset can be used to estimate tree %, but what about all the other LC classes such as sparse vegetation – how did you assign fractions to each PFT? Did you use the ESA CCI cross-walking table for that? And therefore did you consider uncertainties in those very subjective cross-walking values, as described in Hartley et al., 2017.\nResponse: We appreciate the detailed feedback and agree that quantifying uncertainty in RS PFT distributions is a challenging task, as demonstrated by Hartley et al. (2017) and the early comment from the reviewer. Given these complexities, we simplified our approach by relying on reported uncertainties from the literature, rather than developing a framework to recalculate them, which would merit another study. Please also kindly refer to our response to R3C2.2, in which we introduce our updated method to quantify both own and inter-product uncertainty in PFT maps.\nWe note there is a big gap in quantifying the own uncertainty of each RS PFT map, and some maps did not provide any sources of uncertainties in their product at all (Table S3). While MODIS and ESA CCI uncertainties have been reported, existing literature did not consider all uncertainty sources or only for some PFT classes (Hartley et al., 2017; Sulla-Menashe et al., 2019). For MODIS, only the uncertainty of mapping land cover (LC) was provided, lacking the uncertainty for translating from LC to PFT. For ESA CCI, the reported uncertainty accounts for both mapping and translation, but only tailored to major vegetation classes — Tree, Shrub, Grass, and Bare. In our original version, we assigned the uncertainty from the \"Tree\" class to tree-specific PFTs (EBF, ENF, DBF, and DNF), from \"Shrub\" to Shrub PFTs, and similarly, the Grass uncertainties were uniformly applied to cold, cool, and warm grasses.\nAfter some consideration, we think that our original approach may have underestimated RS PFT uncertainties, in agreement with the reviewers’ comments. In this round of revision, we adopted five widely-used RS PFT maps and considered inter-product uncertainty to avoid over-constrained GPP resulting from insufficient consideration of uncertainty when using an individual RS PFT map. We have updated all figures and discussions using the inter-product uncertainty, please see our response to R3C2.\nR3C6: Results – Relationship between PFT GPP and area\nI understand perhaps better than I did in the original review what point the authors are trying to make here. Still, I would say that for some PFTs the linear relationship is not as strong (lower R2 values) and the fact that the R2 values are often lower when they are not forcing the regression line through zero also gives me pause. The logic for forcing the regression through zero is sound, but the authors are also just trying to see what the sensitivity of the model GPP is to the PFT area. In cases where the R2 is much lower, then actually what the non zero constrained regression is telling you is that actually the total GPP is not that sensitive to the area of the PFT as it is for other PFTs. This is especially the case for C3 grasses which are one of the most abundant PFTs worldwide (if you add up the area of the individual C3 categories in this study). So when the authors say in their response to R1C2 that “it is very surprising to see ALL models show EBF has similar productivity per area (i.e., the slopes of the total GPP – area relationships for each PFT in Fig. 1)” I would say I agree, but I also think that’s not the case for all PFTs.\nI think this analysis is interesting though because the differences across PFTs in terms of the strengths of the linear regression (unconstrained R2 values) provide insights into which PFTs have very similar representation across different models (and therefore where the prescribed PFT area is causing most of the difference in total area) versus the PFTs where area plays a role, but clearly other aspects of modeling that PFT are contributing to the spread in total GPP (which then gives lower R2 values).\nSome of the PFTs that clearly have a strong linear relationship between total GPP and area are not that surprising – i.e. C4 grasses have had little attention paid to them compared to C3 grasses, and so it is perhaps predictable that all models simulate a similar GPP per area as their underlying representation is very similar. In contrast, more work has been done on C3 photosynthesis – for example, different models for calculating C3 photosynthesis, different parameter values (some of which have been calibrated against more numerous datasets for C3 grasses), including optimality-based approaches, and including temperature acclimation, etc – and therefore I assume there are differences across models, and that explains the lower R2 values. (I also note that Figure S1 is not cited in the main text). Deciduous needleleaved forests are another example that have had relatively little attention paid to them compared to other PFTs and therefore the R2 values are quite high.\nResponse: Thank you for the comments, which further expanded the comment R3C2.1. We have therefore added a paragraph to specifically discuss the difference in strength of the relationships between PFTs, in L303-307.\n“This framework also helped explain the variability in the strength of GPP-area relationships among PFTs (i.e., R2 and slope values in Fig. 1) – the relatively lower R2 values, such as those for C grasses and crops, may indicate greater inter-model difference in LAI or leaf 3 photosynthetic rates compared to other PFTs.”\nR3C7: Results – Updated land carbon cycle estimation using remotely sensed PFT area\nGiven the issues I mentioned above I am afraid I do not trust the narrow range of uncertainties used for either MODIS or ESA CCI. See my comments in the general section on suggestions for addressing the key caveats above.\nFig. 2: Caption include the fact that the DGVMs are from TRENDY and should say which version, even though this is in the methods (as many people won’t read the methods). In fact I would make it clear what is going into each of these estimates in the figure caption. Which products are going into MLEC for example, just FLUXCOM or a range of different upscaled EC products?\nResponse: Please kindly refer to our answer to R3C2.2 about the uncertainty of RS PFT maps and associated reduction in GPP estimates.\nWe have also clarified the version of TRENDY DGVMs (v9), as well as the version names for ML products in both figure captions and Method section. Specifically, ML includes the EC EC\nRS_METEO output from FLUXCOM – ensemble estimates upscaled using different machine learning methods and meteorological forcings (Jung et al., 2020), as well as X-BASE product from the FLUXCOM-X (Nelson et al., 2024).\nR3C8: Results – Contribution of PFT change to global GPP change\nThis section is really nice and interesting and yet it is a little lost in the overall message about constrained PFT area leading to constrained GPP based on the strength of the linear relationship between PFT GPP and area. I would have liked to see more here, although I appreciate there’s not much space in a letters journal. What is very curious about this is that we expect to see C4 grass outcompete at higher CO2 concentrations, if I’m remembering correctly. This is not directly discussed although the authors mention the potential water savings under elevated CO2. I would suggest the authors expand this though to suggest future studies check that that is indeed what is happening (water savings) and the incrased C3 is not because of incorrect process representation in the models.\nResponse: Thank you. We have expanded the discussion (L189-190) to suggest that “the expansion of C was likely due to its anatomical traits, which confer higher water-use 4 efficiency”, supported by theoretical understanding (Way et al., 2014) as well as global isotope observations (Eggels et al., 2021). The C and C WUE are very different in theory, and C is 3 4 4 expected to be more WUE. If models use conventional Farquhar approach for C and Collatz 3 approach for C , they all should be able to get higher WUE from C than C , though the 4 4 3 magnitude may vary depending on model parameterization.\nR3C9: Discussion\nLines 247-249: “To improve biogeography process in DGVMs, we suggest prioritizing efforts to clear hotspots of biogeography uncertainty (e.g., those identified in Fig. 3) and to develop a more reliable global PFT map” à Yes absolutely! Actually most of these models have no dynamic vegetation nor vegetation demography that is causing any kind of PFT change. They all have prescribed maps. So yes, a more reliable global PFT map is needed for models that do that! And yet, as the authors have shown, the two most widely available maps have very different areas, which, as the authors have also shown, feeds into differences in total GPP (though I just wouldn’t place as much emphasis on the narrowing of uncertainty in that total GPP estimate as the authors do). This is the key message of the paper here I think, but I believe it can be framed more clearly.\nResponse: Following the suggestion from the reviewer, we have avoided putting too much emphasis on narrowing of the GPP uncertainty. Instead, we updated Fig. 3 using five widely-used RS PFT maps and highlighted the uncertainty hotspots in PFT maps. We have rephrased the statements to\n“Our study provided a framework to quantify the magnitude of land carbon cycle using RS PFT maps, however, we must emphasize the considerable difference among RS PFT maps (Fig. 1; Fig. 3). A long line of research has indicated large disagreement between RS PFT maps, and those differences are still challenging to reconcile. The key difference could be relevant to the fraction of specific PFT types in the world (e.g., especially the EBF and C crops). For 3 example, MODIS suggests the EBF area is 21.56 ± 1.23 × 106 km2, but ESA CCI suggests it is only 14.23 ± 2.85 × 106 km2. The discrepancy in RS PFTs can propagate into the global GPP estimates, leading to an approximately 20 Pg C difference (Fig. S2).”\nWe have also included a paragraph to discuss the uncertainty hotspot of RS PFT maps (L212-218).\n“Moreover, we found that RS PFT maps also show large differences in their estimates of PFT total area and spatial distributions (Fig. 3; Table S3). For example, there were large spreads in the estimates of global EBF area (i.e., from 13.97 to 21.56 × 106 km2) and C3C area (from 7.66 to 16.34 × 106 km2) across the five remote sensing products. The hotspots of biogeography disagreement among RS PFT maps include EBF in the eastern U.S. and western Europe, ENF in northern Canada and eastern Russia, DNF in North America and Europe, and SHR and C3G across the globe.”\nR3C10: Line 249-250: I guess I would say be careful with the use of “alternative” in this statement: “Using a consistent biogeography baseline, such as a PFT map from remote sensing as a diagnostic alternative” because these PFT maps are being used in models. I know ORCHIDEE and SDGVM use the same ESA CCI map (which also makes me wonder about the circularity in this approach – I would like the authors to comment on that).\nResponse: We have rephrased the statement to “Using a consistent biogeography baseline in DGVMs can be a first step to reduce PFT-induced uncertainty in land carbon cycle simulations.”\nRegarding the issue of circularity arising from the use of the ESA CCI map, we would like to clarify that although some DGVMs use the ESA CCI map, the implementation of this map varies largely among models:\n1) We note clear differences in the use of ESA CCI product in ORCHIDEE and SDGVM. ORCHIDEE distinguishes among 15 PFTs, whereas SDGVM differentiates only 10. For instance, ORCHIDEE divides C grasses into temperate, tropical, and boreal categories, while 3 SDGVM aggregates all C grasses into a single category. These differences imply distinct 3 PFT-specific parameterizations and structural frameworks, potentially resulting in different GPP rates per unit area.\n2) Even when models employ ESA CCI maps, they often integrate with additional land-cover datasets. SDGVM, for example, merges ESA CCI 2014 land-cover data (Poulter et al., 2015) with LUH2v2h cropland and pasture datasets.\n3) In addition, SDGVM enabled a dynamic vegetation component, while ORCHIDEE used a prescribed static PFT map. The two models also simulated on different spatial resolution – SDGVM at 0.5° and ORCHIDEE at 0.25°.\nGiven the varying approaches for implementation of ESA CCI data across DGVMs, the PFT distributions in these models are not likely to be the same with each other (also shown in our Fig. 1). Additionally, in this round of revision, we have updated the main result based on five RS PFT maps, rather than a single one, to further avoid the issue of overly relying on one RS PFT map.\nR3C11: Lines 256-261: This is only true for dynamic vegetation models, of which there are only 4 in the TRENDY v9 ensemble. Otherwise, photosynthesis rates etc do not feedback on PFT area. So I would revise this section accordingly.\nResponse: Please kindly see our response to R3C3.\nR3C12: Lines 269- : I like the additional emulator analysis and I really appreciate the authors attempting to address this comment from the first round of review. I have one minor comment and one concern:\nMy concern is that I don’t think the larger spread (smaller reduction in uncertainty in EM(DGVMS,MODIS) is because the models have tuned their structure to the wrong location. Or at least, I can’t see how that would be the case. Let’s say C4 grass is in the completely wrong location in many models because the bioclimatic limit used to separate out C3 and C4 grasses is wrong. If so, what does that imply, that their Vcmax parameter for example is adjusted to the climate of the wrong location to match in situ flux tower data? That seems convoluted to me, although I think maybe I see the argument. I’d be willing to bet though that most PFTs are not in the completely wrong location (mostly based on prescribed PFT maps). It would mean major biases in parameters or model structure to account for that. What I think is maybe more likely, or explains more of the difference between EM(DGVMs,MODIS) and DGVM_MOD is what I’ve said throughout this review – that the uncertainty reduction from constraining with MODIS (or ESA CCI) PFT map area is too high. I would guess that if you had more realistic uncertainties on your PFT maps, your constrined GPP would also be about 40-50% as in EM(DGVMs,MODIS). It still shows a big importance of PFT map, but again, I think it just cannot explain the majority of the spread across models. Still, this nice emulator result shows that differences simply due to thye PFT map are a greater issue than anyone would bet.\nMy minor comment is that I am not sure I would describe it as a difference between biogeography and biogeochemistry (Responses to R1C7 and R2C3 and new text describing Fig. S5) (with my issue being mostly with the use of biogeochemistry). If the location of the PFTs is wrong (e.g., Bias type 1 in Response to R2C3), then that’s also a biogeography issue right? Or am I misunderstanding? If the location is wrong then it could be that the bioclimatic envelopes for example are not right, and I would call that a biogeography “process” or limitation (that is intimately tied to area). On the other hand, you may have an issue with temperature optima, which would prohibit certain PFTs from growing in an area. I would call that a biogeochemistry (or plant physiology actually) process (could also be related to phenology thresholds or ecohydrology if not enough water was available). So it’s just a question of wording but I perhaps would more simply say that what you said in response to R2C3 – that there are two types of biases related to area and location, and that biases related to area only were accounted for by the analysis in Figure 1. But Figure 3 showed the second bias related to location distribution, and so you have done this analysis. Basically, I think the use of biogeochemistry here is confusing and also misleading given that for most models any “biogeochemistry/plant physiology” etc process is not feeding back on the distribution anyway. If the distirbution is wrong for most models then it’s because the prescribed PFT map is wrong (except for the dynamic vegetation models).\nResponse: We appreciate the discussion from the reviewers on the interpretation of our Fig. S5 (the emulator analysis, now Fig. S7), and how much the figure tells us about the relative importance of biogeography and biogeochemistry in the uncertainty of modelling GPP.\nFigure S7. The spread of global GPP estimates from original DGVMs and DGVMs forced by MODIS PFT map. The three group of models are 1) DGVMs – the original group of DGVMs, with biased PFT area and location; 2) EM(DGVMs, MODIS) – the emulators of DGVMs forced with MODIS PFT maps, meaning the PFT area and locations are both “corrected”; 3) DGVM RS – the constrained DGVMs simulations, meaning only the PFT area was “corrected” but location remained biased.\nPerhaps it is good for us to establish some baseline understanding first. We quote reviewer’s comment below.\n1. “I’d be willing to bet though that most PFTs are not in the completely wrong location (mostly based on prescribed PFT maps).” – this is a complex case. Below we show the maps of the disagreement of PFT maps across DGVMs (Fig. S6), and they can be really different and thus wrong (at least for some)! For example, a few models predict abundant EBF distributions (> 10%) in parts of Siberia and ENF in tropical South America and Africa.\nFigure S6. The difference in the distribution of PFTs in DGVMs, quantified by another method. Here we presented the number of DGVMs in each grid cell for each PFT. The number of DGVMs where PFT percentage larger than 10% in the same grid, with higher values indicating a greater consistency. The abbreviations for PFTs correspond to those in Figure 1.\n2. “that the uncertainty reduction from constraining with MODIS (or ESA CCI) PFT map area is too high”. Originally, we reported 86% reduction, and now it is 75% after we used five RS PFT maps to constrain GPP. Even with this updated result, biogeography still emerges as a main factor controlling the uncertainty of the land carbon cycle modelling.\n3. Thus, the intriguing finding illustrated by Fig. S5 (now Fig. S7) remains valid. When we correct area (DGVM ), uncertainty in GPP simulations reduced by 75%; when we correct RS both area and location (EM(DGVMs, MODIS)), uncertainty reduced by 46%. It means that, if we place PFT at the right location, it may not immediately reduce the uncertainty, because the new parameters coming with the replaced PFT may actually inflate the uncertainty. We do feel without a proper MIP and scenario run, we are unable to provide a clear picture, but our (best) guess is that the parameters of the correct PFT for that location might have some issues.\nIn this analysis, since we are discussing inter-model differences using multi-year average, the issue whether parameters feedback on PFT distribution is not a direct concern, given their static PFT maps incurred a large difference in GPP estimates already. Please also see our reply to R3C3.\nR3C13: Reply to response to R1C3 in response to reviewers: I would perhaps disagree that PFTs can be directly linked to spectral signatures. In fact, it is very hard to distinguish trees from grasses, even in hyperspectral data and even when you add LiDAR data (as my PhD student Rubaya Pervin found in drylands given there’s also a strong spectral signature from dry grasses and/or biocrusts). In the past this has prompted the suggestion of using “optical functional types” (Ustin and Gamon, 2010) where global vegetation would be classified based on the clustering of spectral signatures. To my knowledge, no one has really tried this at regional or global scale. I would bet this would look different to the PFT maps we currently use and would be hard to match with anything we represent in models. So the question remains, how can we use satellite data to infer plant functional type fractional cover? The fusion of tree cover fraction in the ESA CCI product is a great advance. But it also depends on what is counted as a tree (height thresholds etc).\nResponse: Thank you for sharing your thoughts on remote sensing-based PFT mapping. We apologize for the somewhat broad statement on “spectral signature”, and we have updated it to “spectral signature and its variation” (such as phenology, snow cover) to reflect what is really being implemented in many RS PFT mapping work.\nWe encourage community input to improve PFT mapping, though addressing this grand challenge is beyond the scope of our study here – we intend to explore the global vegetation biogeography in the future, not only based on remote sensing, but also based on ecological theories and the insights they might provide to PFT distributions.\nR3C14: Reply to response to R1C6 in response to reviewers: Yes this is very interesting to see the range of estimates, and also the reporting of more recent studies that all have higher estimates. In fact, the higher mean value of GPP as constrained with MODIS is one of the most interesting things about this paper (although this is somewhat dampened by the fact that ESA CCI maintains the typical value of around 120PgC), although other methods have a higher value as the authors show in Fig. 2a. As I am signing this review I apologise for this shameless self-citation, but a paper in which I constrained GPP using SIF data (one of the earlier, simpler approach using just one model) also found a posterior (post parameter optimization) global value of just over 160PgC (MacBean et al., 2018). Please do not feel you need to include this is a reference, it is more for your interest. Since I wrote that paper in the back of my mind has been the thought that yes ok we constrained the high positive bias in ORCHIDEE (pre optimization it was above 190PgC), but we have not constrained it enough. I hadn’t quite processed until I read the authors response to the original review that there are indeed multiple studies pushing at this higher value (and hadn’t seen the Liu paper). As the authors say, the value of 120-130PgC has been talked about for Beer’s 2010 Science paper if not before and is the number that is most often quoted. It has been 130PgC in the the GCB for at least the last two years. So I think it’s great that the authors have moved away from a “correct” magnitude and instead show the uncertainties. Still, when you put all these studies together with the new study presented here, it does feel like the evidence is mounting for a higher value.\nResponse: We thank the reviewer for sharing her earlier study, and we are glad to see another piece of evidence supporting a higher GPP value. We cited it to further support our argument that based on the evidence so far, it is likely global GPP is greater than we thought before. The jewel is out there again. Our study may motivate the community to think about why we should have a higher GPP than we thought before.\nMinor points:\nR3C15: Line 29: would suggest “biogeography-climate” instead of “biogeography – climate” interactions. First time I read it I thought the words after the hyphen were making a different point.\nResponse: Thanks. We have revised this accordingly.\nR3C16: Lines 54-56: another limitation for MODIS GPP (or any GPP derived from the production efficiency/LUE method) and the machine learning approaches I believe is the lack of accounting for the effect of photosynthesis on LUE (De Kauwe et al., 2016).\nResponse: Yes, general LUE-based models, such as the one used in MODIS GPP product, do not account for elevated CO₂ (De Kauwe et al., 2016) and their considerations of climate impacts are quite empirical (Medlyn, 2011). We have now added this point in L59.\nR3C17: Line 71: would maybe not use the word “prescribe” here (to be pedantic) – maybe just that in process-based models these processes are specific to each PFT ?\nResponse: We have removed the word ‘prescribe’ and updated the sentence to “the current framework for improving land carbon cycle simulations in DGVMs implies the importance of biogeography, given that changes in model parameters, structures, and even the emerging climate sensitivities of simulated carbon fluxes can also be driven by PFT change, as DGVMs largely regard these processes as PFT-specific”.\nR3C18: I cannot see CLM5 in Fig. 1j? I’m guessing that’s why you have axes upper limits that are much higher than the main cluster of models\nResponse: CLM5 does not appear in Fig. 1j because, although the TRENDY v9 dataset indicates that CLM5 provides C crop data, upon careful examination of the output files we 3 found the data for C crop to be entirely NaN across all pixels. Therefore, we were unable to 3 include CLM5 C crop in this analysis. 3\nR3C19: I like the new additional analysis presented in Fig. S6. Nice!\nResponse: We are pleased to hear that the additional analysis is helpful.\nR3C20: Data Availability Statement needs to include where you got the TRENDY v9 GPP and PFT from as it’s actually not on the global carbon budget page that links to Mike O’Sullivans Github (https://mdosullivan.github.io/GCB/). Also, would you consider providing the emulators in a publicly available repository for people to test other PFT maps?\nResponse: We obtained the TRENDY v9 data from the project coordinator, Dr. Stephen Sitch and have now revised this statement to,\n“DGVM simulations from TRENDY v9 were obtained from Dr. S. Sitch (s.a.sitch@exeter.ac.uk) and are now available at https://mdosullivan.github.io/GCB/.”\nMoreover, we have now made the emulator code and associated models publicly available via the Zenodo repository at https://doi.org/10.5281/zenodo.15752958.\nReferences\nDe Kauwe, M. G., Keenan, T. F., Medlyn, B. E., Prentice, I. C. and Terrer. C. (2016) Satellite based estimates underestimate the effect of CO2 fertilisation on net primary productivity. Nature Climate Change, 6, 892-893.\nMacBean, N., Maignan, F., Bacour, C., Lewis, P., Peylin, P., Guanter, L., ... & Disney, M. (2018). Strong constraint on modelled global carbon uptake using solar-induced chlorophyll fluorescence data. Scientific reports, 8(1), 1973.\nUstin, S. L., & Gamon, J. A. (2010). Remote sensing of plant functional types. New Phytologist, 186(4), 795-816.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We greatly appreciate your review and have endeavoured to use all comments to improve our study. Please kindly find our detailed point-by-point responses to reviewers 3 and 4 altogether above.\nReferences:\nBuchhorn, M. et al., 2020. Copernicus global land service: land cover 100m: version 3 globe 2015-2019: algorithm theoretical basis document. Zenodo.\nDantas De Paula, M. et al., 2025. Including the phosphorus cycle into the LPJ-GUESS dynamic global vegetation model (v4.1, r10994) – global patterns and temporal trends of N and P primary production limitation. Geoscientific Model Development, 18(7): 2249-2274.\nDe Kauwe, M.G., Keenan, T.F., Medlyn, B.E., Prentice, I.C. and Terrer, C., 2016. Satellite based estimates underestimate the effect of CO2 fertilization on net primary productivity. Nature Climate Change, 6(10): 892-893.\nEggels, S., Blankenagel, S., Schön, C. and Avramova, V., 2021. The carbon isotopic signature of C4 crops and its applicability in breeding for climate resilience. Theoretical and Applied Genetics, 134(6): 1663-1675.\nFatichi, S., Pappas, C., Zscheischler, J. and Leuzinger, S., 2019. Modelling carbon sources and sinks in terrestrial vegetation. New Phytologist, 221(2): 652-668.\nFriedlingstein, P. et al., 2023. Global Carbon Budget 2023. Earth System Science Data, 15(12): 5301-5369.\nHarper, K.L. et al., 2023. A 29-year time series of annual 300 m resolution plant-functional-type maps for climate models. Earth System Science Data, 15(3): 1465-1499.\nHartley, A.J., MacBean, N., Georgievski, G. and Bontemps, S., 2017. Uncertainty in plant functional type distributions and its impact on land surface models. Remote Sensing of Environment, 203: 71-89.\nHu, Z. et al., 2018. Joint structural and physiological control on the interannual variation in productivity in a temperate grassland: A data-model comparison. Global Change Biology, 24(7): 2965-2979.\nJung, M. et al., 2020. Scaling carbon fluxes from eddy covariance sites to globe: synthesis and evaluation of the FLUXCOM approach. Biogeosciences, 17(5): 1343-1365.\nKeenan, T.F. et al., 2023. A constraint on historic growth in global photosynthesis due to rising CO2. Nature Climate Change, 13(12): 1376-1381.\nLiu, H. et al., 2020. Annual dynamics of global land cover and its long-term changes from 1982 to 2015. Earth Syst. Sci. Data, 12(2): 1217-1243.\nMedlyn, B.E., 2011. Comment on “Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 Through 2009”. Science, 333(6046): 1093-1093.\nNelson, J.A. et al., 2024. X-BASE: the first terrestrial carbon and water flux products from an extended data-driven scaling framework, FLUXCOM-X. Biogeosciences, 21(22): 5079-5115.\nPoulter, B. et al., 2015. Plant functional type classification for earth system models: results from the European Space Agency's Land Cover Climate Change Initiative. Geoscientific Model Development, 8(7): 2315-2328.\nPrentice, I.C. et al., 1992. Special Paper: A Global Biome Model Based on Plant Physiology and Dominance, Soil Properties and Climate. Journal of Biogeography, 19(2): 117.\nRogers, A. et al., 2017. A roadmap for improving the representation of photosynthesis in Earth system models. New Phytologist, 213(1): 22-42.\nSulla-Menashe, D. and Friedl, M.A., 2018. User guide to collection 6 MODIS land cover (MCD12Q1 and MCD12C1) product. Usgs: Reston, Va, Usa: 1-18.\nSulla-Menashe, D., Gray, J.M., Abercrombie, S.P. and Friedl, M.A., 2019. Hierarchical mapping of annual global land cover 2001 to present: The MODIS Collection 6 Land Cover product. Remote Sensing of Environment, 222: 183-194.\nWang, X. et al., 2021. Beneficial effects of nitrogen deposition on carbon and nitrogen accumulation in grasses over other species in Inner Mongolian grasslands. Global Ecology and Conservation, 26: e01507.\nWay, D.A., Katul, G.G., Manzoni, S. and Vico, G., 2014. Increasing water use efficiency along the C3 to C4 evolutionary pathway: a stomatal optimization perspective. Journal of Experimental Botany, 65(13): 3683-3693.\nZhang, X. et al., 2024. GLC_FCS30D: the first global 30 m land-cover dynamics monitoring product with a fine classification system for the period from 1985 to 2022 generated using dense-time-series Landsat imagery and the continuous change-detection method. Earth Syst. Sci. Data, 16(3): 1353-1381.\nZhu, J. et al., 2025. Changing patterns of global nitrogen deposition driven by socio-economic development. Nature Communications, 16(1): 46."}]

Decision letter