Representation of the terrestrial carbon cycle in CMIP6

Simulation of the carbon cycle in climate models is important due to its impact on climate change, but many weaknesses in its reproduction were found in previous models. Improvements in the representation of the land carbon cycle in Earth system models (ESMs) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) include the interactive treatment of both the carbon and nitrogen cycles, improved photosynthesis, and soil hydrology. To assess the impact of these model developments on aspects of the global carbon cycle, the Earth System Model Evaluation Tool (ESMValTool) is expanded to compare CO₂-concentration- and CO₂-emission-driven historical simulations from CMIP5 and CMIP6 to observational data sets. A particular focus is on the differences in models with and without an interactive terrestrial nitrogen cycle. Overestimations of photosynthesis (gross primary productivity (GPP)) in CMIP5 were largely resolved in CMIP6 for participating models with an interactive nitrogen cycle but remaining for models without one. This points to the importance of including nutrient limitation. Simulating the leaf area index (LAI) remains challenging, with a large model spread in both CMIP5 and CMIP6. In ESMs, the global mean land carbon uptake (net biome productivity (NBP)) is well reproduced in the CMIP5 and CMIP6 multi-model means. However, this is the result of an underestimation of NBP in the Northern Hemisphere, which is compensated by an overestimation in the Southern Hemisphere and the tropics. Carbon stocks remain a large uncertainty in the models. While vegetation carbon content is slightly better represented in CMIP6, the inter-model range of soil carbon content remains the same between CMIP5 and CMIP6. Overall, a slight improvement in the simulation of land carbon cycle parameters is found in CMIP6 compared to CMIP5, but with many biases remaining, further improvements of models in particular for LAI and NBP is required. Models from modeling groups participating in both CMIP phases generally perform similarly or better in their CMIP6 compared to their CMIP5 models. This improvement is not as significant in the multi-model means due to more new models in CMIP6, especially those using older versions of the Community Land Model (CLM). Emission-driven simulations perform just as well as the concentration-driven models, despite the added process realism. Due to this, we recommend that ESMs in future Coupled Model Intercomparison Project (CMIP) phases perform emission-driven simulations as the standard so that climate–carbon cycle feedbacks are fully active. The inclusion of the nitrogen limitation led to a large improvement in photosynthesis compared to models not including this process, suggesting the need to view the nitrogen cycle as a necessary part of all future carbon cycle models. Possible benefits when including further limiting nutrients such as phosphorus should also be considered.

1 Introduction

Earth system models (ESMs) simulate the climate system by interactively coupling physical general circulation models of the atmosphere, ocean, land, and cryosphere with biogeochemical and biophysical cycles (Jones, 2020). The Coupled Model Intercomparison Project (CMIP; Meehl et al., 2000) was established to facilitate a consistent comparison between different ESMs through the use of common forcings and a uniform output structure in order to better understand past, present, and future climate. The newest phase, CMIP6 (Eyring et al., 2016a), provides a large ensemble of model simulations and includes 23 CMIP6-Endorsed Model Intercomparison Projects (MIPs) which facilitate a better analysis of specific scientific questions. Every new phase of CMIP sees additional models and improved model components, making a comparison to previous phases vital to determine if known systematic biases have been reduced and model weaknesses were identified and overcome (Eyring et al., 2019) or if increased model realism through the inclusion of additional processes introduces new biases. Increasing the process realism of models, for example, by replacing time-invariant observation-based fields with interactive prognostic ones, while having a neutral impact on the present-day performance of the model, can be viewed as a successful step in model improvement. This is particularly important for more self-consistent climate projections. Bock et al. (2020) assessed atmospheric variables of the CMIP6 ensemble and compared them to CMIP3 and CMIP5 output. They find that for temperature, precipitation, water vapor, and zonal wind speed, many long-standing biases remain in the multi-model mean, but individual models and high-resolution versions of models show some improvements for temperature and precipitation. In this study, we expand this assessment to the carbon cycle.

The carbon cycle is an important part of ESMs due to the key role of anthropogenic emissions of carbon dioxide (CO₂) in driving climate change (IPCC, 2021). The land carbon cycle dominates the uncertainty in the global carbon cycle (Canadell et al., 2021). The Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) was largely supported by model simulations from the CMIP5 ensemble, and CMIP6 models were an important input to AR6 (IPCC, 2021). It is important to evaluate these models and document their changes compared to CMIP5. The focus of this study is on the results of the CMIP6 historical simulations that are split into simulations with prescribed greenhouse gas concentrations and those with prescribed CO₂ emissions. Simulations driven by carbon emissions require an interactive carbon cycle to determine the distribution of natural and anthropogenic carbon fluxes across the land, marine, and atmospheric reservoirs instead of relying on prescribed atmospheric CO₂ concentrations (Friedlingstein et al., 2014). Therefore, only emission-driven simulations have fully active climate–carbon cycle feedbacks vital for self-consistent future projections. Sanderson et al. (2024) thus petition for the use of emission-driven simulations to be prioritized in CMIP7. This leads to the further importance of the evaluation of the carbon cycle as more models will be required to implement and improve their interactive carbon cycle. An analysis of idealized 1 % CO₂ increase per year (1pctCO₂) simulations for CMIP5 and CMIP6 models with respect to carbon concentration and carbon–climate feedbacks has been carried out by Arora et al. (2020). They found that while these feedback parameters have not changed significantly between CMIP5 and CMIP6, the land feedback parameters are weaker for models including a nitrogen cycle coupled to the carbon cycle. Davies-Barnard et al. (2020) documented the development and inclusion of modeling of the terrestrial nitrogen cycle in ESM land surface schemes and showed how this affects the response to elevated CO₂ across models. They find that these models show a more accurate response in the tropics than high-latitudes compared with observed responses. Gier et al. (2020) investigated the atmospheric CO₂ concentrations for emission-driven CMIP5 and CMIP6 models and found that while CMIP6 models show an improvement in reproducing the observations compared to CMIP5, during the period of the satellite observations (2003–2014) the growth rate is overestimated, and the seasonal cycle amplitude is underestimated in both CMIP5 and CMIP6 relative to observations. Furthermore, the model spread in simulated atmospheric CO₂ (∼ 45 ppmv in CMIP5, ∼ 35 ppmv in CMIP6) was found to remain many times larger than the observational uncertainty of under 1 ppmv (parts per million by volume) over this period.

Anav et al. (2013a) investigated the land and ocean carbon cycle for CMIP5 historical model simulations. While most models reproduced the main climatic variables and their seasonal evolution correctly, weaknesses were found in their ability to reproduce more specific biogeochemical fields such as a general overestimation of photosynthesis. Consequently, these were some of the main areas tackled in updating and enhancing the land carbon models of the ESMs for CMIP6, including the addition of a coupled nitrogen cycle and nitrogen limitation, which can limit the rates of carbon cycling through vegetation and soil, as well as photosynthesis updates and soil hydrology improvements (e.g., Danabasoglu et al., 2020; Delire et al., 2020; Wiltshire et al., 2021). Jones et al. (2023) drew on an expert assessment of regional carbon budgets to evaluate the terrestrial carbon cycle in CMIP6 models. They found that the multi-model mean performs well in most regions for most variables (both carbon fluxes and stocks), but individual models have strengths and weaknesses. This paper assesses the impacts of the model improvements and additions in CMIP6, especially the impact of an additional coupled nitrogen cycle, and whether they help overcome some of the weaknesses identified in CMIP5. Expanding on the analysis from Anav et al. (2013a) for CMIP5, this study uses CMIP6 concentration- and emission-driven historical simulations to compare to the land carbon cycle in CMIP5 and identify possible improvements in the newer model generation. In this study, we restrict ESM evaluation to data sets with global extent.

Section 2 describes the data used in this study, while Sect. 3 contains the analysis of different carbon cycle variables in the comparison between CMIP ensembles and observations. Long-term trends and the seasonal cycle of carbon cycle variables are considered, and the analysis is concluded with a performance metrics plot evaluating the climatological seasonal cycle of the models with different observational data sets. Section 4 summarizes the results and conclusions.

Law et al. (2017)Ziehn et al.20172020Ziehn et al. (2019a, b)Swart et al. (2019a)Swart et al. (2019e, b)Swart et al. (2019a)Swart et al. (2019d, c)Danabasoglu et al.2020Danabasoglu (2019b)Danabasoglu et al.2020Danabasoglu (2019a)Cherchi et al. (2019)Lovato and Peano2020Lovato et al. (2022)Lovato et al. (2021)Séférian et al. (2019)Séférian (2018, 2019)Döscher et al. (2022)EC-Earth Consortium2021b2021aDöscher et al. (2022)EC-Earth Consortium2019Dunne et al. (2020)Krasting et al.2018a2018bVolodin et al. (2018)Volodin et al. (2019a)Volodin et al. (2017a, b)Volodin et al. (2019b)Boucher et al. (2020)Boucher et al. (2018)Hajima et al. (2020b)Hajima et al.20192020aNeubauer et al.2019aTegen et al.2019Neubauer et al. (2019b)Mauritsen et al. (2019)Wieners et al.2019a2019bYukimoto et al. (2019a)Yukimoto et al.2019c2019bSeland et al. (2020)Seland et al. (2019b, a)Seland et al. (2020)Bentsen et al. (2019)Sellar et al. (2019)Tang et al. (2019b, a)Park et al. (2019)Park and Shin (2019)Lee et al. (2020)Lee and Liang (2020)

Table 1 CMIP6 models analyzed in this study. In the comments column, D stands for models including dynamic vegetation and N for models including nitrogen cycles. Models for which emission-driven simulations are also analyzed are marked in bold.

Ji et al. (2014)Arora et al. (2011)Hurrell et al. (2013)Dunne et al. (2012, 2013)Dunne et al. (2012, 2013)Collins et al. (2011)Team et al. (2011)Collins et al. (2011)Team et al. (2011)Volodin et al. (2010)Bao et al. (2012)Qiao et al. (2013)Dufresne et al. (2013)Dufresne et al. (2013)Watanabe et al. (2011)Watanabe et al. (2011)Giorgetta et al. (2013)Giorgetta et al. (2013)Yukimoto et al. (2011)Tjiputra et al. (2013)

Table 2 CMIP5 models used in this study, with notations as in Table 1.

2 Data and tools

2.1 CMIP simulations

Model simulations from both CMIP Phase 6 (Eyring et al., 2016a) and Phase 5 (Taylor et al., 2012) are used, with Tables 1 and 2 listing model characteristics such as their atmosphere and land model components, in addition to their main references. A more comprehensive summary of the land model components of the CMIP models is given in Appendix A. Models were selected due to their availability on Earth System Grid Federation (ESGF) nodes for the considered variables (see the “Code and data availability” section).

This study focuses on historical simulations which aim to reproduce the observed climate since the pre-industrial times. They span 1850–2005 (CMIP5) and 1850–2014 (CMIP6). Both experiments with prescribed greenhouse gas concentrations (concentration driven; historical) and prescribed CO₂ emissions (emission driven; esmHistorical and esm-hist for CMIP5 and CMIP6, respectively) are considered but evaluated separately and compared to each other. Models participating in the emission-driven simulations, marked in bold in the tables, use their interactive carbon cycle to determine the distribution of natural and anthropogenic carbon fluxes across the land, marine, and atmospheric reservoirs instead of relying on prescribed atmospheric CO₂ concentrations (Friedlingstein et al., 2014).

Very few CMIP5 models had a coupled carbon–nitrogen cycle. While the BNU-ESM model included carbon–nitrogen interactions, they were turned off for the CMIP5 model simulations as the nitrogen cycle had not been fully evaluated (Ji et al., 2014). Therefore, a nitrogen cycle was included in 2 out of 18 CMIP5 models (CESM1-BGC and NorESM1-ME) which both use the CLM4 land model and in 15 out of 23 CMIP6 models spread over six different land models, with the Community Land Model (CLM) in different versions accounting for 8 CMIP6 models (v5: CESM2, CESM2-WACCM, NorESM2-LM, and NorESM2-MM; v4.5: CMCC-CM2-SR5 and CMCC-ESM2; v4: SAM0-UNICON and TaiESM1). The other land models in CMIP6 with a coupled nitrogen cycle are LPJ-GUESS (EC-Earth3-CC and EC-Earth3-Veg), JSBACH (MPI-ESM-1-2-HAM and MPI-ESM1-2-LR), CABLE+CASA-CNP (ACCESS-ESM1-5), JULES-ES (UKESM1-0-LL), and Visit-e (MIROC-ES2L). This shows a large bias towards the CLM land model in CMIP6 which needs to be considered while analyzing the multi-model mean (MMM).

To facilitate a direct comparison of CMIP5 and CMIP6 data in figures containing temporal means, only data up to 2005 representing the end of the CMIP5 historical simulations are considered. Unless stated otherwise, figures use mean data over the time period 1986–2005. Only one realization per model is used, as different ensemble members perform similarly to each other with respect to the carbon cycle (Gier et al., 2020), and using an ensemble mean would lead to an under representation of the internal variability present in individual ensemble members (Anav et al., 2013a). Multi-model means (MMMs) were computed separately for each CMIP phase and experiment combination, as well as an additional distinction between models with and without interactive nitrogen models, and are computed on the monthly gridded data for which models are regridded to a common 2°×2° grid. MMMs are not weighted according to the interdependence of the models and model components or according to their performance relative to observational products.

While we split models into groups only dependent on the presence of an interactive nitrogen cycle in this study, vegetation dynamics is another important process for ESM comparison. Models interactively simulating vegetation cover may simulate trees or grasses in the wrong areas compared to models using observational land cover maps, thus impacting variables with a strong relation to land cover, such as leaf area index (LAI) or gross primary productivity (GPP). While models with prescribed land cover may show better LAI in the present day, they cannot predict future changes in vegetation cover or their impact on regional climate and carbon processes. For reference, Tables 1 and 2 note models with dynamic vegetation with a D in the comment column.

2.2 Reference data

A large range of observations and reanalysis data sets have been used to assess model performance. These data sets are listed in Table 3, along with their main reference(s), their source, the variables used, and their temporal coverage. Both observational and reanalysis data sets will be referred to as “observations” from here on, in contrast with the results from the CMIP model simulations. The longest observational records are derived from reanalyses, while satellite observations only provide data from the late 20th century onward. Since most reference data sets do not come with observational uncertainty, a common approach is to use several reference data sets per variable where available, as noted in Seiler et al. (2022). This approach is also taken in this study.

Rödenbeck (2005)Chevallier et al. (2005, 2010)Chevallier (2013)Friedlingstein et al. (2022)Global Carbon Project (2021)Jung et al. (2019)Jung et al. (2011)Yuan et al. (2007)Liang et al. (2021)Zhu et al. (2013)Cao et al. (2023)Gibbs (2006)Wieder (2014)Hugelius et al. (2013)

Table 3 Reference data sets used in this study. Data sets in bold are the main reference data set, and those in italics are the alternate references for Figs. 15–17.

For the leaf area index (LAI), we use the LAI3g product (Zhu et al., 2013) that provides global monthly gridded data starting in the year 1981. It has been generated using an artificial neural network based on data from the Advanced Very High Resolution Radiometer (AVHRR) and the Moderate Resolution Imaging Spectroradiometer (MODIS). Furthermore, we also use the newly released GIMMS LAI4g data set (Cao et al., 2023), which is based on the same satellite data as LAI3g but employs a newer Normalized Differential Vegetation Index (NDVI) data set base which removes the effects of the satellite orbital drift and AVHRR sensor degradation which plagued many other LAI data sets. Furthermore, LAI4g uses a large number of high-quality Landsat LAI samples to increase the spatiotemporal consistency of the data set. Last, the Global Land Surface Satellite (GLASS; Liang et al., 2021) is a product suite with 12 products and from which we employ both the leaf area index (LAI) and gross primary productivity (GPP) products. Similar to LAI3g and LAI4g, GLASS is based on data obtained from AVHRR and MODIS. As newer GLASS data products only use MODIS and thus start from 2000, this paper uses a previous GLASS version (v4.0) which includes AVHRR data and thus starts in 1981 for LAI and 1982 for GPP, respectively. GLASS LAI uses general regression neural networks trained on preprocessed reflectance data of an entire year to estimate the 1-year LAI profile for each pixel. The LAI product is one of the variables used to estimate GLASS GPP with an eddy covariance light use efficiency model. Both GLASS products are available on a 0.05° grid with a frequency of 8 d.

Another GPP product, MTE (Jung et al., 2011), provides global monthly gridded data starting in 1982. It uses an upscaling of data from the FLUXNET eddy covariance tower network based on the model tree ensemble (MTE) approach. Similarly, the FLUXCOM product (Jung et al., 2019) is also based on an upscaling of FLUXNET site level observations but additionally incorporates a larger variety of machine learning methods and also includes remote sensing (from MODIS) and meteorological data. Here, we use a global monthly gridded version of FLUXCOM (starting in 1980) from the RS+METEO setup. Due to the assumption of an unchanging average CO₂ level, both MTE and FLUXCOM data are known to have an unrealistic non-existent trend (0.01 PgC yr⁻² globally) (Anav et al., 2015). Thus, trend analysis on GPP should exclude these data sets.

The main data set for the land–atmosphere carbon flux (net biome productivity (NBP)) is the Jena CarboScope (version sEXTocNEET_v2020) product (Rödenbeck, 2005), which provides global daily gridded data starting from the year 1957. This data set provides surface–atmosphere CO₂ fluxes based on atmospheric measurements calculated from an atmospheric transport inversion. The inversion used here (NEE-T inversion) involves a regression of interannual net ecosystem exchange (NEE) anomalies against air temperature anomalies (T). In total, Jena CarboScope uses data from 156 atmospheric measurement sites distributed across the entire globe. The alternative data set for the land–atmosphere carbon flux is a further inversion product from the Copernicus Atmosphere Monitoring Service (CAMS; Chevallier et al., 2005, 2010; Chevallier, 2013). CAMS provides global gridded data on a monthly resolution starting in 1979 (other temporal resolutions are also available). The inversion product we use here (v20r2) is based on surface measurements from more than 100 sites. A third data set used for comparing the global annual mean NBP is the Global Carbon Project (GCP; Friedlingstein et al., 2022), which estimates the global carbon budget using several observations and models. It provides estimates for emissions from fossil fuel combustion and industrial processes, emissions from land use change, the atmospheric CO₂ growth rate, ocean sink, land sink, cement carbonation sink, and the budget imbalance from combining all these terms. The land–atmosphere carbon flux for GCP has to be calculated by subtracting the land use change emissions from the residual land sink. The land sink was obtained from averaging the results from 17 dynamic global vegetation models (DGVMs) which reproduce the observed mean total land uptake of the 1990s and is given with an uncertainty of ± 0.5 PgC yr⁻¹ on average. The land use change emissions are estimated from the average of three bookkeeping models with an uncertainty of ± 0.7 PgC yr⁻¹, making it one of the only data sets with direct estimations for uncertainties.

For each of the remaining carbon cycle variables, only a single reference data set is taken into account. For the vegetation and soil carbon pools, the NDP-017b (Gibbs, 2006) and HWSD + NCSCD (Wieder, 2014; Hugelius et al., 2013) products are used, respectively. Both data sets provide global gridded annual data for the single year of 2000. NDP-017b uses an updated database that extends the methodology of Olson et al. (1985), who developed a global carbon stock map of above- and belowground biomass using 20 years of field investigations, consultations, and literature analysis and the more up-to-date land cover conditions of the Global Land Cover Database (GLC2000). The Harmonized World Soil Database (HWSD) uses large volumes of regional and national soil information to create an empirical data set that provides soil parameter estimates for topsoil (0–30 cm) and subsoil (30–100 cm). Similar to Varney et al. (2022), we combine the HWSD data set with the Northern Circumpolar Soil Carbon Database (NCSCD; Hugelius et al., 2013) to complement the HWSD data in the polar region. It uses data on soil order coverage to calculate soil organic carbon content and mass with 1778 pedon data. Wherever overlap between the two data sets occurs, the NCSCD data are chosen.

2.3 ESMValTool

The analysis in this paper was produced using the Earth System Model Evaluation Tool (ESMValTool) version 2 (Righi et al., 2020; Eyring et al., 2020; Lauer et al., 2020; Weigel et al., 2021; Schlund et al., 2023). ESMValTool is an open-source community diagnostics and performance metrics tool which has been developed to routinely evaluate ESMs contributing to CMIP and compare them with other ESMs, predecessor versions, as well as observations. Since its first release (Eyring et al., 2016b), ESMValTool has been updated for increased performance in its core functionality to deal with the increased data volume of CMIP6 and now features full traceability and reproducibility through provenance, as well as new and updated diagnostics and metrics, which can be applied to many models and variables. Available diagnostics cover a large range of scientific topics and are described in three papers. They include large-scale diagnostics for quasi-operational and comprehensive evaluation of ESMs (Eyring et al., 2020), diagnostics for extreme events, regional model and impact evaluation and analysis (Weigel et al., 2021), and diagnostics for emergent constraints and analysis of future projections (Lauer et al., 2020). A new suite of recipes has been developed to cover the work of this study, as well as some improvements on previous diagnostics for the carbon cycle available in ESMValTool. This facilitates the evaluation of the carbon cycle in future studies, including the analysis of upcoming CMIP7 simulations that can be easily compared to CMIP5 and CMIP6 to assess improvements.

3 CMIP model performance

General climate variables, such as temperature and precipitation have a large influence on the carbon cycle. It is therefore important to assess how well these variables are simulated by the ESMs. If they are well reproduced but carbon cycle variables are not, then it is likely due to a poor representation of processes specific to the carbon cycle, while a poor performance in the physical variables makes an attribution of the cause of poor performance in the carbon cycle variables more difficult. The CMIP5 and CMIP6 ensembles have been assessed compared to observations by Flato et al. (2013) and Eyring et al. (2021), respectively. A detailed analysis was also done by Bock et al. (2020), and references therein, who compare the surface temperature, pressure, precipitation, radiation, and clouds of CMIP3, CMIP5, and CMIP6 historical simulations for annual means. The CMIP6 models show better correlations for these variables than the CMIP5 models for all parameters, with smaller improvements for variables such as temperature which were already well represented in previous CMIP phases. However, the model spread is not significantly reduced but instead largely remains the same. Here, we expand the analysis to the carbon cycle. However, it should be mentioned that many carbon cycle processes are affected by physical variables on much smaller timescales, such as the timing of precipitation throughout the day or if surface temperatures falling below the freezing point at any time of the day may impact the growth of plants more than their monthly means suggest. This study uses monthly mean data, which does not resolve many of these important events and thus does not investigate the impact of physical variables further, as any analysis would still be lacking many possible impacts of sub-frequency effects. Future studies using higher-frequency data should put more emphasis also on the physical drivers of carbon processes, as well as the role of extreme events.

3.1 Leaf area index

The leaf area index (LAI) is the ratio of one-sided leaf area per unit of ground area (Anav et al., 2013a) as a measure of the canopy structure. Models use LAI to calculate the photosynthetic uptake of the total canopy, also known as gross primary productivity (GPP). While LAI is an important building block for the carbon cycle, it was also one of the weaknesses of the carbon cycle in the CMIP5 ensemble and tended to be overestimated (Anav et al., 2013a, b).

Carbon uptake by land follows a pronounced seasonal cycle, with CO₂ removed from the atmosphere through plant photosynthesis and released back through plant and soil respiration. With LAI describing the canopy structure and more plants thriving in summer, it is strongly linked to the seasonal cycle of atmospheric CO₂. The seasonal cycles for LAI for CMIP5 and CMIP6 MMMs for concentration- and emission-driven simulations are shown in Fig. 1 and split into models with (Ncycle) and without interactive nitrogen cycle (non-Ncycle), as well as different regional means, namely global, Northern Hemisphere (20–90° N; NH), Southern Hemisphere (20–90° S; SH), and tropics (20° N–20° S). From here on, concentration-driven simulations will be denoted by c, such as CMIP5c and CMIP6c, and emission-driven simulations by e (CMIP5e and CMIP6e). The BNU-ESM and MRI-ESM1 CMIP5 models were removed from the MMM due to featuring an unrealistically high mean LAI, almost doubling the LAI of the reference data and the other models in the SH and the tropics. A common mask is applied to all data sets, which includes all missing values in any data set to allow for direct comparison between the models and reference data sets. This increases the LAI compared to unmasked regional means, as missing values are more common in desert and mountainous regions with low LAI (Fig. 3). The chosen reference data sets LAI3g, LAI4g, and GLASS agree well across all regions, which was to be expected, as they are all based on the same raw satellite data from AVHRR and MODIS. Xiao et al. (2017) found that the GLASS product outperformed other products, which included LAI3g, when compared to LAI from high-resolution reference maps. As such, GLASS, which is also the reference data set with the largest coverage, will be considered the main reference data set for our analysis of LAI.

Figure 1 Seasonal cycle of leaf area index (LAI) for a climatological mean of 1986–2005 for different regional averages, namely global, Northern Hemisphere (20–90° N), Southern Hemisphere (20–90° S), and tropics (20° S–20° N). The reference data sets (LAI3g, solid line; LAI4g, dashed line; GLASS, dotted line) are shown in black, while the MMMs for the different phase–experiment combinations are denoted by blue for CMIP5 and red for CMIP6, with darker colors for the emission-driven simulations (dark blue for CMIP5e; dark red for CMIP6e) and lighter colors for concentration-driven simulations (light blue for CMIP5c; light red for CMIP6c). MMMs derived from models with coupled nitrogen cycle (Ncycle) are dashed, while solid lines represent MMMs of models without a coupled nitrogen cycle (non-Ncycle). The shading represents the standard deviation of the MMMs, with vertical hatching for models without and horizontal hatching for models with coupled nitrogen cycle. For comparison with the reference data, which contain many missing values, a common mask was applied to all data sets, removing values where any data set is missing a value.

There is a strong seasonal cycle in the NH, dominated by high-latitude vegetation in Eurasia and North America. The NH seasonal cycle is the dominant contribution to the global mean due to the higher relative land fraction in the higher latitudes of the Global North compared to the Global South. All models overestimate LAI in all regions, but the CMIP6 models reproduce the reference data better than the CMIP5 models. According to Anav et al. (2013b), the overestimation in the mid-latitudes is likely partly due to a wet bias and its control on soil moisture, a saturation of satellite instrumentation, and missing parameterizations of disturbances. The annual mean precipitation wet bias is minimally reduced in CMIP6 (Bock et al., 2020), and new parameterizations such as nutrient limitations through interactive nitrogen cycle have been introduced in some models, leading to a reduced LAI in the CMIP6 MMMs compared to CMIP5. The seasonal cycle in the NH is reproduced, but while the CMIP5 non-Ncycle models reproduce the amplitude well with a positive offset of approximately 0.7 m² m⁻², the CMIP6 non-Ncycle models are better at reproducing the peak value. Both CMIP5 Ncycle models (CESM1-BGC and NorESM1-ME) use the CLM4 land model with known issues regarding LAI, such as underestimating LAI in dry regions due to elevated CO₂ and overestimating LAI in moist regions (Lee et al., 2013), as well as an unrealistically strong nitrogen limitation (Wieder et al., 2019), hindering plant growth. This leads to these models showing a larger overestimation in LAI in both the Southern Hemisphere and tropic regions dominated by moist rainforests. Additionally, the seasonal cycle amplitude is strongly reduced in the NH, while the mean LAI is larger than for the reference data sets. Both CMIP phases show a weakness in simulating the end of the growing season, shown by the later decline in LAI in winter, which also leads to a smaller seasonal cycle amplitude, consistent with the findings of Park and Jeong (2021). The drawdown in autumn signifying the end of the growing season is smaller in Ncycle models compared to non-Ncycle models. The differences between the concentration- and emission-driven simulations are small, with models participating in both simulations having very similar results (individual models not shown). Larger differences occurring here and in later analysis between the concentration- and emission-driven simulations are likely due to the different subset of models in the historical simulations and not due to the experiment design.

Figure 2 shows the mean and trend of LAI averaged over 1986–2005 and depicts phase–experiment simulations with one type of marker each for a better overview. Exact numbers for all data sets are found in the Supplement (Tables S1–2). The errors given in the tables refer to the standard deviation of the mean as a measure for the interannual variability (IAV), while the standard error in the trend is the error in the linear regression calculating the trend. For the individual reference data, LAI3g and LAI4g agree well in mean and trend, while GLASS agrees with their mean but shows a significantly higher trend in all regions, leading to the large trend error bar. Xiao et al. (2017) analyzed the trend of several LAI products for different biome types for 1982–2011 and found GLASS to have significantly higher trends in savannahs and shrubs compared to LAI3g but lower trends in deciduous broadleaf forests and evergreen needleleaf forests, and even a negative trend for deciduous needleleaf forests, while grasses, cereal crops, and evergreen broadleaf forest trends are similar for GLASS and LAI3g. This results in larger trend differences in the SH (0.009 m² m⁻² yr⁻¹) and tropics (0.010 m² m⁻² yr⁻¹) than the NH (0.003 m² m⁻² yr⁻¹). The MMM global means show a significant overestimation compared to the reference data, with the overestimation of the CMIP5 MMMs of 0.7 m² m⁻² reduced by half for CMIP6 MMMs, as previously seen in Fig. 1. This large improvement for the CMIP6 models is not related to the difference between Ncycle and non-Ncycle models, as their LAI MMM global means are comparable for CMIP6. Due to the large difference in trends in the reference data sets, the global mean trends of all CMIP MMMs lie within the range of the reference data. The model trends range between slightly negative and strong positive. Unlike the mean, the LAI trend does not show a strong difference between CMIP phases or any other grouping method we employed. The CMIP6 models only show a slightly smaller range in trend compared to the CMIP5 models, but more CMIP5 models have a lower trend compared to the reference data sets than CMIP6 models.

Figure 2 Mean and trend of LAI computed over 1986–2005 for different regions, namely global, Northern Hemisphere (20–90° N), Southern Hemisphere (20–90° S), and tropics (20° S–20° N). The mean of the reference data sets (LAI3g, LAI4g, and GLASS) is denoted by a black star (“”), with error bars for the standard deviation. Models for phase–experiment combinations are shown with a single symbol each, with the blue smaller than (“<”) sign for CMIP5 and the red greater than (“>”) sign for CMIP6, while darker colors show emission-driven simulations, and lighter colors show concentration-driven simulations, as well as Ncycle models being denoted with a filled symbol. MMMs are depicted with crosses (“x”) in the color assigned to their phase–experiment combination.

Figure 3 Global maps of LAI averaged over 1986–2005 for all reference data sets, as well as the mean and range between lowest and highest values per grid cell of the reference data sets.

The other regions reflect these global MMM comparisons as well. CMIP6 means are closer to the reference data than CMIP5 in all regions, although still overestimating LAI, and agree well with each other no matter the experiment or Ncycle status. The only exception is the NH CMIP6c Ncycle MMM, which shows a larger mean than the other three CMIP6 groupings. This is due to the CMCC-ESM2 and CMCC-CM2-SR5 models, which show a much higher LAI in the NH compared to the reference data but fit well in the other regions. These two models use the CLM4.5 land model, which Li et al. (2022) found to have a far longer peak growing season and to overestimate LAI in boreal forests compared to MODIS reference data, consistent with our results. The CLM5 models (CESM2, CESM2-WACCM, NorESM2-LM, and NorESM2-MM) perform much better in the mean LAI than CLM4.5 in the NH but are still overestimating LAI compared to the reference data. The trends in the NH for the reference data are on the lower end of a large model range. Only the CMIP5 non-Ncycle MMMs agree with the reference data, while the other MMMs show a larger trend than the reference data but are comparable to each other. CMIP6 MMM LAI mean and trend values agree well with the reference data in the tropics and the SH, while the CMIP5 MMMs overestimate the mean by ≈ 0.7 m² m⁻² in the SH and ≈ 0.6 m² m⁻² in the tropics for non-Ncycle, as well as ≈ 1.2 m² m⁻² in the SH and ≈ 1.1 m² m⁻² in the tropics for Ncycle MMMs, but show a similar trend to LAI3g, LAI4g, and CMIP6 MMMs. The larger mean LAI in the CLM4 (CMIP5 Ncycle MMMs) can be traced back to the overestimation of LAI in the moist regions mentioned before. Some models show a significant negative trend in LAI in the SH and the tropics, resulting in a globally negative trend even with a positive trend in the NH.

Maps of the LAI reference data are found in Fig. 3 but without the common mask to see the different coverages. Coverage of the different reference data sets varies a lot due to different quality control criteria and algorithms, with most missing values found in desert or mountainous regions such as the Sahara and the Himalayas. Additionally, a mean of the reference data and the range of the reference data per grid cell is shown, along with the global mean of the values in the upper-right corner. GLASS has a larger coverage over desert and mountainous regions, which are regions with low plant coverage and thus low LAI, resulting in a lower global mean LAI of 1.36 m² m⁻² compared to 1.67 and 1.71 m² m⁻² from LAI3g and LAI4g, respectively. This underlines the importance of the common mask used for Figs. 1 and 2 to obtain comparable results. The different data sets show the same pattern of the LAI distribution, with the absolute values ranging between 1 and 6 m² m⁻², while the differences are below 2 m² m⁻², with the largest difference occurring in tropical rain forests and northern high latitudes – the regions with the largest absolute LAI values. For a grid cell bias comparison of the different model groupings (Fig. 4), a combined reference data set was computed as the mean of the other reference data sets. Due to the different coverages, some areas are only calculated from the GLASS data, while others are an average of all three data sets. The range of values per grid cell going into the combined data set is plotted in the lower right of Fig. 3. The largest difference occurs in the areas with larger LAI, i.e., the tropical rainforests and then boreal forests, with a global mean average range of 0.38 m² m⁻². For the MMM bias maps shown in Fig. 4, hatching is added where the MMMs agree with the reference mean within the MMM standard deviation. CMIP5 Ncycle MMMs show the issue of CLM4 in overestimating LAI in wet regions, with LAI in tropical rainforests almost doubling the reference value, while drier regions show a significant negative bias. While this results in a global mean bias of 0.66 to 0.67 m² m⁻², smaller than 0.89 m² m⁻² for CMIP5e non-Ncycle MMMs, it is still the worst-performing model grouping when considering a grid cell basis due to its extreme biases in both directions. The hatching showing the agreement can be ignored in this case, as only two models contributed to the MMM standard deviation. The CMIP5e non-Ncycle MMM shows a strong overestimation across the northern latitudes besides Greenland. This can be attributed to the GFDL-ESM2G and GFDL-ESM2M models, which are known to have established coniferous trees in areas which should contain tundra or cold deciduous trees, as in its vegetation spin-up only coniferous trees are allowed to grow in cold regions, but grasses or deciduous trees which would have a lower LAI are not (Anav et al., 2013b). While the GFDL models show this problem in both CMIP5c and CMIP5e, due to the larger number of MMMs contributing to the concentration-driven simulations, their effect is reduced. In CMIP6, GFDL-ESM4 still has a large positive LAI bias throughout these areas, but it is significantly reduced compared to its CMIP5 predecessor. The second prominent overestimation is around the tropical rainforests, where models like BNU-ESM and MRI-ESM1, and to a lesser extent, the GFDL models also extend the LAI hotspot to larger areas around it compared to the reference mean. In CMIP5c, HadGEM2-CC and HadGEM2-ES also show this overestimation. CMIP5c non-Ncycle MMM shows no special bias patterns, but instead a general overestimation in almost all areas with hatching is present throughout the globe. CMIP6 MMMs show a reduced mean bias of less than half the CMIP5 overestimation, with almost no pronounced patterns and a bias reduction in all areas, with the largest improvements found in the northern high latitudes. The largest bias is in Southeast Asia for CMIP6c Ncycle MMMs, which can be tracked to the CMCC-ESM2 and CMCC-CM2-SR5 models, and this makes the mean bias of the CMIP6c Ncycle MMM higher than that of CMIP6c non-Ncycle MMM. Otherwise, the bias pattern looks similar for CMIP6 Ncycle and non-Ncycle models.

While CMIP6 LAI has improved compared to CMIP5, especially a significantly reduced mean bias, a general overestimation of LAI remains, along with issues of correctly reproducing the length of the growing season in the Northern Hemisphere and a large model spread in mean and trend LAI. Neither the introduction of an interactive nitrogen cycle nor the comparison between emission- and concentration-driven simulations shows large differences in the overall quality of the CMIP6 simulations for LAI.

Figure 4 Global maps of LAI bias for 1986–2005 with respect to the reference data set mean shown in Fig. 3. The panels show the MMMs of the models with (left) and without (right) coupled nitrogen cycle for the different phase-experiment combinations. The hatching represents the areas where the MMM of the models and reference mean agree within the MMM standard deviation.

3.2 Gross primary productivity

Gross primary productivity (GPP) represents the CO₂ uptake on land due to photosynthesis. This was one of the biggest weaknesses of the CMIP5 ensemble, with most models overestimating photosynthesis, as well as leaf area index (Anav et al., 2013a). The seasonal cycle of GPP (Fig. 5) shows good agreement between the GLASS and MTE reference data, while the FLUXCOM data show a lower GPP in all regions, as well as a shorter growing season in the Northern Hemisphere. All models reproduce a similar shape of the NH seasonal cycle to the GLASS and MTE data in both model generations and experiments. As found in Anav et al. (2013a), the CMIP5 non-Ncycle models overestimate GPP in all regions, while the CMIP5 Ncycle models strongly underestimate the peak of the seasonal cycle in the NH. CMIP6 models perform better than CMIP5, but while the CMIP6 non-Ncycle models still overestimate the GPP peak in summer similarly to CMIP5 for the NH, the Ncycle models show a very good agreement for both CMIP6c and CMIP6e. Nitrogen limitation is stronger in northern latitudes through boreal forests and tundra (Du et al., 2020) compared to tropical and subtropical forests, which are more limited by phosphorus. However, from the CMIP models used in this study, only ACCESS-ESM1-5 includes an interactive phosphorus cycle. It therefore makes sense that Ncycle models show a decreased GPP compared to non-Ncycle models in the NH and are closer to the reference data. The model spread remains large in CMIP6, despite being smaller for Ncycle models, which is denoted by horizontal hatching. Following LAI, there is no strong and discernible seasonal cycle in the SH or the tropics, but the CMIP6 models are closer to the mean GPP than the CMIP5 models, with lower values for Ncycle models.

Figure 5 As in Fig. 1 but for gross primary production using GLASS, FLUXCOM, and MTE reference data. Additionally, the regional GPP is calculated as the area-weighted sum instead of the mean of the grid cells used for LAI.

The temporal mean and linear trend of the spatial sums for GPP during the time period 1986–2005 is shown in Fig. 6. MTE and FLUXCOM data are known to have an unrealistic non-existent trend (0.01 PgC yr⁻² globally) due to the assumption of an unchanging average CO₂ level (Anav et al., 2015) in these data sets. As such, the model trend should not be compared to the trend of these two reference data sets, and we have omitted these data sets from the calculation of the reference trend. The mean and trend for GPP of all regions of these data sets, along with all other numerical values from the plot, can be found in Tables S3–S4. The trend of GLASS is closely linked to the high trend of LAI GLASS, which is one of the main influences on GPP. The reference data sets agree well with respect to the mean GPP, with the largest difference being a lower mean for FLUXCOM in the NH. The models show a large range, with MRI-ESM1 CMIP5e as an outlier, which shows an even larger mean GPP in all regions. The CMIP5 models are on the higher side of this range, with the Ncycle CMIP5 models much lower due to their underestimation in the NH as seen in Fig. 5. The CMIP6 Ncycle MMMs agree very well with the reference data, while the CMIP6 non-Ncycle MMMs show a larger mean GPP. The global trend of GLASS is positive, which is well matched by the non-Ncycle CMIP6c MMM with the other MMMs showing a smaller trend. In the NH, more CMIP5 models overestimate the mean GPP than CMIP6 models. The MMMs span a large range, with the Ncycle MMMs showing a lower mean GPP than the non-Ncycle MMMs, which are overestimating GPP compared to the reference data. While the trends of the models are centered around the GLASS trend, there is a large model range with outliers for MRI-ESM2-0 in both its CMIP6c and CMIP6e runs. The CMIP6 non-Ncycle MMM shows a higher trend than the other MMMs due to the MRI-ESM2-0 outliers. In the SH mean GPP, the CMIP6 MMMs match the reference data well, with lower values for Ncycle than for non-Ncycle MMMs. The CMIP5 MMMs are slightly above these, with the non-Ncycle CMIP5e having a much larger value due to the outliers of MRI-ESM1 and FIO-ESM, which both have values well above 100 PgC yr⁻¹. The MMMs underestimate the trend compared to the reference data. The distribution in the tropics is very similar to the SH. The reference data trend is underestimated by the MMMs. In summary, the Ncycle MMM shows a better performance than the non-Ncycle MMM in the NH, while it shows a slight underestimation in the tropics and a similar performance in the SH. The model spread over the trend in GPP stays similar throughout the model generations, with the mean trend being largely consistent with the GLASS reference data set in all regions but underestimated everywhere but in the NH.

Figure 6 As in Fig. 2 but for gross primary production using GLASS, FLUXCOM, and MTE reference data.

As for LAI, the coverage of the GLASS data is larger than for the other reference data, with missing values for FLUXCOM and MTE mainly found over the Sahara and the Himalayas (Fig. 7). GLASS shows a larger GPP in the tropical rainforests and boreal forests, explaining the larger mean GPP seen before. The reference mean data set has hardly any missing values left, and the largest differences in the data sets are at places with the highest GPP, with discrepancies in the areas bordering the hotspots of the rainforests and boreal forests. Even though the GPP bias maps (Fig. 8) for Ncycle CMIP6e MMMs have a global mean bias almost 2 orders of magnitude smaller than the CMIP6e non-Ncycle MMMs, they show the same pattern of overestimation in wet regions and underestimation in dry regions found for LAI (Fig. 4), underlining the strong influence of LAI on GPP. The CMIP5 non-Ncycle MMMs also shows similar patterns to the LAI bias maps, but the overestimation in the areas around the tropical rainforests is strongly reinforced, while the previous strong overestimation of LAI in the northern high latitudes for the non-Ncycle CMIP5e MMM is reduced. The CMIP5c non-Ncycle MMM additionally shows a strong underestimation at the northeastern coast of South America. The global mean bias for CMIP5 non-Ncycle MMMs lies at 0.36 kgC m⁻² yr⁻¹ for CMIP5e and is reduced by half for CMIP5c. This bias is further reduced to approximately 0.1 kgC m⁻² yr⁻¹ for CMIP6 non-Ncycle MMMs, which show similar bias patterns to CMIP5 non-Ncycle MMMs, but there is an overall reduction in the bias patterns, with a larger reduction in the savannah regions of Africa. The global mean bias for the CMIP6 Ncycle MMMs is further reduced to 0.001 kgC m⁻² yr⁻¹ for CMIP6e and a negative bias of −0.04 kgC m⁻² yr⁻¹ for CMIP6c. Both show a reduction in the Northern Hemisphere bias of the open shrublands, turning some into a negative bias, as well as southwest Africa, while the slight overestimation in North America remains, as well as the underestimation at the northeastern part of South America. This is summarized in Fig. 9, which shows the zonal sums of the reference data and the MMMs. Unlike the seasonal cycle and scatterplots shown before, a common mask is not applied here, but instead values are masked out if more than 15 % of the pixels over land in a latitude band are missing for a given data set. The large overestimation of non-Ncycle CMIP5e can be seen, which is reduced in non-Ncycle CMIP5c, with both cycles showing a peak slightly north of the Equator which is not seen in the reference data and which is due to the overestimation of the shrublands south of the Sahara. The CMIP6 non-Ncycle MMMs show a much better approximation across all latitudes, with a slight reduction in the bias in the NH, but they still show a significant overestimation in the tropics. This is remedied in the CMIP6 Ncycle models, which show a very good agreement with the reference data across all latitudes and now with slight underestimations at high latitudes.

Figure 7 Global maps of GPP averaged over 1986–2005 for all reference data sets (GLASS, FLUXCOM, and MTE), as well as the mean and range between lowest and highest values per grid cell of the reference data sets. The number in the top right denotes the global GPP flux.

Figure 8 Global maps of GPP bias for 1986–2005 with respect to the reference data set mean shown in Fig. 7. The panels show the MMMs of the models with (left) and without (right) coupled nitrogen cycle for the different phase–experiment combinations. The hatching represents the areas where the MMM of the models and reference mean agree within the MMM standard deviation, while the number in the top right denotes the global mean bias.

3.3 Land–atmosphere flux

The net carbon flux from the atmosphere into the land (net biome productivity – NBP) characterizes the balance between carbon uptake due to photosynthesis and carbon release by respiration, as well as other processes like fires and de- and afforestation. Positive values of NBP denote carbon uptake by land. The CMIP6 EC-Earth models (EC-Earth3-CC and EC-Earth3-Veg) are excluded from the MMM for NBP because they show a very strong land source in December due to the annual update of land use transitions (Döscher et al., 2022), thus influencing the seasonal cycle. For the consistency of models contributing to the MMM, these models are not only excluded from the MMM calculation of the seasonal cycle but all NBP MMMs. However, their individual data points are included in the scatterplot. MIROC-ESM and MIROC-ESM-CHEM are also removed from the MMMs due to rapidly changing fluxes both temporally and spatially. This occurs due to the stochastic dynamic vegetation model (Watanabe et al., 2011) and shows stronger effects in the tropics when combined with the land use change emission module and the high climate sensitivity. While the influences on the regional means or climatologies are smaller than for EC-Earth3, they can be seen in mean map plots very well.

The global seasonal cycle for NBP (Fig. 10) is dominated by the Northern Hemisphere, with almost no discernible cycle in the Southern Hemisphere or the tropics. There is generally a good agreement between the two inversions, but the CAMS inversions show a larger seasonal cycle in the SH and tropics of approximately 6 PgC yr⁻¹, where CarboScope shows no clear seasonal cycle. In the NH, and due to its large contribution to the total also globally, CAMS has a higher NBP at the start of the year and NBP to a lesser degree at the end of the year, where CarboScope shows a larger negative NBP and thus a carbon sink. The models agree with the carbon sink of CarboScope in these months but have a weaker carbon sink (higher NBP) in NH autumn. The CMIP5 Ncycle MMMs have a smaller seasonal cycle amplitude compared to any of the other MMMs and the reference data, which are carried over from GPP. The other MMMs reproduce the seasonal cycle well, while the non-Ncycle CMIP6e MMM is shifted late by a month, showing possible issues with the start and end of the growing season. In the SH and tropics, where CarboScope found no significant cycle and CAMS had a slightly larger one in the tropics, the Ncycle models follow the shape and timing of the CAMS data, while the non-Ncycle models have a seasonal cycle shifted to 2 months earlier. There is no significant difference between CMIP5 and CMIP6 or between the c and e experiments in the MMM.

Figure 9 Area-weighted zonal sums of gross primary production and the reference data sets of GLASS, FLUXCOM, and MTE. No common masking is applied, but latitudes were set to missing if more than 15 % of the land grid cells contained missing data. The hatching depicts the MMM standard deviation, with a horizontal hatching for models with an interactive nitrogen cycle and vertical hatching for models without an interactive nitrogen cycle.

Figure 10 As in Fig. 1 but for the land–atmosphere carbon flux and the reference data sets of CAMS and CarboScope.

The temporal mean and trend of spatially summed NBP are shown in Fig. 11, with specific values given in Tables S5 and S6. CAMS shows a larger mean NBP compared to CarboScope in all regions but the tropics, which is consistent with the NBP averages from Seiler et al. (2022), using this data set, who found its NBP to be larger than comparable data sets and model results. GCP as a global average is only available as a reference data set for the global panel. The models show a far larger range for the global mean NBP than the reference data, with outliers for the INM-CM CMIP6c models showing a far larger mean NBP than the other models. The Ncycle CMIP5 means have a negative mean NBP, while the other MMMs show a better agreement with the reference data. The Ncycle MMMs have significantly smaller mean NBPs than the non-Ncycle MMMs. The relatively good overall agreement of the models' global mean NBP with the reference data does not hold for the different regions. Most models and all MMMs simulate a lower carbon sink in the Northern Hemisphere when compared to the inversions, with Ncycle models generally showing a smaller mean NBP but no large discernible differences between the different groupings. Conversely, while the inversions estimate both the Southern Hemisphere and the tropics to be a slight carbon source due to deforestation, the MMMs, with the exception of the CMIP5 Ncycle, show a carbon uptake by land in these regions. The large values for the non-Ncycle CMIP6c MMM are again due to the overestimation of the INM-CM4-8 and INM-CM5-0 models, but as their mean NBP in the NH is not a large outlier, we did not remove these from the MMM. The underestimation in the NH combined with the overestimation in the SH and tropics leads to the good global agreement of the total carbon sink. This is in agreement with the findings from IPCC AR6 (Eyring et al., 2021; Canadell et al., 2021). The inclusion of a nitrogen cycle and therefore the inclusion of nitrogen limitations on CO₂ fertilization was expected to address this discrepancy of the distribution of the carbon sinks (Canadell et al., 2021), but the data do not support this, as the Ncycle MMMs do not show a different performance to the non-Ncycle MMMs in CMIP6. While the models show a large range of trends, the MMMs agree well with the reference data, and this continues in the other regions as well. While N limitation is not expected to be substantial for the present day, it represents a major limitation on future land–carbon uptake (Zaehle et al., 2015), and thus, its inclusion a major advance in being able to robustly simulate future carbon balance of the terrestrial carbon cycle.

Figure 11 As in Fig. 2 but for the land–atmosphere carbon flux and the reference data sets CAMS and CarboScope. The GCP data are a globally averaged time series and thus only appear in the global plot.

As seen before, the reference data sets show different means, trends, and slightly different seasonal cycles in the different regions. For a more detailed look, Fig. 12 shows maps of the reference data. In this case, the hatched area is the area where the data sets agree on the sign or within a margin of half the bin size of the contour plot. While in large parts of the globe the data sets agree in sign, there are significant differences. In North America, CAMS shows a much larger carbon sink than CarboScope, which instead shows some carbon sources along the west coast and throughout South America. In CAMS, South America is split into a much stronger carbon source in the Amazonian Rainforest and a carbon sink in the southern part. The data sets also disagree in Europe, which is a carbon source according to CAMS but a sink in CarboScope. The literature found Europe to be a carbon sink for the first part of the 21st century (Ciais et al., 2013; Reuter et al., 2014), using different reference data sets. This would support the CarboScope data set, but due to the different time frames considered, it is not definitive. CAMS also sees a carbon sink in tropical Africa, where CarboScope has a slight carbon source. In Southeast Asia, CarboScope has a large carbon sink, where CAMS shows a neutral carbon flux. The area of most agreement is in the NH as a large carbon sink, where there is no deforestation and a bit of afforestation. The Amazon region was a large carbon sink, which is becoming a source due to deforestation (Gatti et al., 2021). CAMS sees the Amazon as a strong carbon source, while CarboScope shows a smaller source with a sink in the northwestern region. Keenan and Williams (2018) found inverse models to show South America as a carbon source and a carbon sink, and it is thus a hotly debated area. Kou-Giesbrecht et al. (2023) attribute the weak agreement between CarboScope and CAMS to differences in the inversion models and atmospheric CO₂ measurements used, with larger differences at latitudes with smaller land areas. Due to the difference between the observational data sets, maps depicting the bias with respect to the reference mean shown for the other variables so far have been omitted. Instead, the area-weighted zonal sums are plotted in Fig. 13 for a comparison of the MMMs with both reference data sets. The reference data sets disagree for almost all latitudes, thus making the model comparison to the reference data in these regions not very meaningful. The area where the reference data is in most agreement is in the northern high latitudes (50–80° N), where both of them show a strong carbon sink, about double that shown in both CMIP5 and CMIP6 models. While the issues in CLM4 and thus the CMIP5 Ncycle models are clearly visible (similarly to LAI and GPP), the other MMMs show similar NBP in all latitudes.

Figure 12 Similar to Fig. 7 but for the land–atmosphere carbon flux. Additionally, the forward-slash hatching symbolizes areas where the reference data sets agree on the sign or where the difference is smaller than the size of one bin of the contour plot.

Figure 13 As in Fig. 9 but for the land–atmosphere carbon flux with CAMS and CarboScope reference data.

3.4 Carbon stocks

Another large uncertainty in CMIP5 was the amount of carbon stored in soil and vegetation. This leads to large uncertainties in land use change emissions, which are important for quantifying cumulative emissions, as well as climate mitigation strategies (Friedlingstein et al., 2023). Varney et al. (2024) investigated the carbon–climate feedbacks of soil and vegetation carbon and found soil carbon to be the dominant response of the land surface, highlighting the need to reduce the uncertainty in carbon storage to better quantify future changes in the climate system. Figure 14 shows a scatterplot of global vegetation carbon against the combined soil and litter carbon. Note that some models (CanESM5-CanOE CMIP6c, GFDL-ESM4 CMIP6e, INM-CM4-8 CMIP6c, INM-CM5-0 CMIP6c, FIO-ESM CMIP5e, CanESM2 CMIP5c, and inmcm4 CMIP5c) did not have data on the ESGF nodes for soil or vegetation carbon and are thus missing from the carbon stock analysis. BNU-ESM CMIP5c and CMIP5e show a far larger vegetation carbon than the other models and are thus removed in the calculation of the mean. Additionally, CLM5 and thus CESM2, CESM2-WACCM, NorESM2-LM, and NorESM2-MM include a full vertical soil profile. For these models, the cSoilAbove1m variable is used for better comparison with the other models, as done in Varney et al. (2022).

Figure 14 Scatterplot of global mean vegetation and soil carbon over 1986–2005 for models, with observations from NDP (vegetation carbon) and HWSD + NCSCD (soil carbon). As in Fig. 2, filled symbols denote models with the nitrogen cycle.

The large spread in the global carbon stocks still remains in CMIP6, as shown in Fig. 14 with values for each data set listed in Table S7. In CMIP5, vegetation carbon was spread between 335 PgC (MPI-ESM-LR CMIP5c) and 802 PgC (GFDL-ESM2M CMIP5e), with the outlier of BNU-ESM even reaching values above 1200 PgC. The spread has only marginally been reduced in CMIP6 from a range of 333 PgC (EC-Earth3-Veg CMIP6c) to 724 PgC (CNRM-ESM2-1 CMIP6e). The reference data are at a value of 478 PgC, in the lower range of the models, with the MMMs ranging between 465 PgC (non-Ncycle CMIP6c) and 547 (Ncycle CMIP5e). The spread in soil carbon is even larger, with a CMIP5 range from 513 PgC (CESM1-BGC CMIP5c) up to 3092 PgC (MPI-ESM-MR CMIP5c). The overestimation by MPI-ESM is due to its decomposition parameterization depending on soil moisture and showing maxima in continental dry lands. In CMIP6 MPI-ESM1-2, the soil carbon model was changed to YASSO, which simulates more plausible soil carbon content (Mauritsen et al., 2019). The spread in soil carbon was not significantly reduced in CMIP6 with a range from 514 PgC (GFDL-ESM4 CMIP6c) to 2913 PgC (CMCC-ESM2 CMIP6c), with the reference value for HWSD + NCSCD at 1561 PgC. The CMIP5 Ncycle models have a soil carbon on the lower end of the range, consistent with the CMIP6 models TaiESM1 and SAM0-UNICON which also use CLM4. The CMIP5 Ncycle MMMs are on the very low end of the range, with 532 and 534 PgC for CMIP5 and CMIP5e, respectively, while the other MMMs are closer to the reference data and range between 1197 PgC for non-Ncycle CMIP6c and 2040 PgC for non-Ncycle CMIP5e. While the CMIP6 Ncycle MMMs are closer to the reference data, no significant improvement due to the inclusion of the interactive nitrogen cycle can be seen when considering the whole spread of the models. This is consistent with Wang et al. (2022), who found that changing models from C to CN coupling often result in lowered ecosystem storage but due to different parameterizations simulate similar carbon pools. Varney et al. (2022) suggest that much of the uncertainty in carbon stocks is due to the simulation of below-ground processes; this is backed up by the differences in soil carbon being much greater than in GPP and thus implicating differences in simulated residence times (Carvalhais et al., 2014; Todd-Brown et al., 2014). For a more in-depth discussion, we would like to refer to dedicated studies, such as Varney et al. (2022) and Wei et al. (2022). Furthermore, Varney et al. (2023) found that while the CMIP6 future soil carbon projections have a lower model spread compared to CMIP5, the structure of soil carbon models within CMIP6 ESMs has likely contributed towards this reduction.

3.5 Overall model performance

Here we assess the overall performance of CMIP5 and CMIP6 models with respect to carbon cycle variables. Figures 15 and 16 show a portrait plot similar to Gleckler et al. (2008) for NBP, GPP, and LAI for all four considered regions (global, Northern Hemisphere, Southern Hemisphere, and the tropics), with the normalized relative root mean square difference in the model simulations to reference data sets. Most models have similar scores when compared to the different observations for GPP and LAI, showing that the inter-model spread in CMIP6 is larger than the observational uncertainty in these variables. For NBP, however, models can have different scores to the considered reference data, which is due to the difference in the reference data found in the previous sections. CMIP6 models on average perform much better than CMIP5 models, with models that had a predecessor in CMIP5 improving on their CMIP5 performance in almost all variables, such as GFDL and IPSL in all variables and CESM and NorESM in LAI, with the exception of CanESM, which shows a reduced performance for NBP in CMIP6. Large improvements can be found in all variables going from CMIP5 to CMIP6, especially in LAI, with the exception of the models using the older CLM4 land component (SAM0-UNICON and TaiESM1) in CMIP6 and GPP, which were previously identified as weaknesses in CMIP5. Only the MRI-ESM2-0 model shows a bad performance in both these variables. As mentioned before, dynamic vegetation in models plays a large role in their ability to simulate variables directly related to it, with models interactively simulating vegetation cover (labeled D in Tables 1 and 2), showing a below-average score for LAI and GPP root mean square difference (RMSD).

Figure 15 Normalized relative space–time root mean square difference (RMSD) of the climatological seasonal cycle with respect to reference observations for CMIP5 (left) and CMIP6 (right) concentration-driven simulations compared to reference data sets. The normalization is done relative to the ensemble median of both CMIP5 and CMIP6 models, with positive values (red) denoting a higher RMSD and thus worse performance, while negative values (blue) denote a lower RMSD than the ensemble median and thus a better performance. The considered time periods depend on the start of the reference data (see Table 3) and end in 2005 to accommodate the end of the CMIP5 data. When using two observational references, a diagonal split is introduced, with the default reference data set being shown on the lower right, while the alternate data set is used for the top-left triangle. The default and alternate reference data sets are marked in Table 3 and are as follows: LAI with LAI4g (main) and GLASS (alt), GPP with FLUXCOM (main) and GLASS(alt), and NBP with CarboScope (main) and CAMS (alt). Models with a nitrogen cycle are marked with blue labels. MMMs for both Ncycle and non-Ncycle models were added, with the models which were excluded from MMMs due to various issues, as stated in the previous sections, also removed from the MMM here. Note that the MMMs were calculated on the climatologies prior to the calculation of the RMSDs, so over- and underestimations can cancel each other out. This is the standard for the performance metric plot implemented in ESMValTool and kept for consistency.

Figure 16 As in Fig. 15 but for emission-driven simulations.

Models which were remarked upon in the previous sections as having good or bad agreement with the observations in specific areas, such as the NBP problems in December for EC-Earth3 have RMSDs that reflect these statements, making this metric a well-suited measure for overall performance. Most models have similar RMSDs in the different regions, with the global value reflecting a mean of the different regions. There does not seem to be a qualitative difference between Ncycle and non-Ncycle models as a whole, but the MMMs perform better than any individual model. The better performance of the MMMs is mathematically expected as long as the assumption that both observations and model simulations draw from the same distribution holds true (Christiansen, 2018). The global NBP, LAI, and GPP are also found in Fig. 42 of Chap. 3 of the IPCC AR6 (Eyring et al., 2021), which showed not only carbon cycle variables but also other land, ocean and atmosphere variables averaged over 1980–1999 for comparison across models from CMIP3 to CMIP6. Their results are compared to reference data from JMA-TRANSCOM for NBP, LAI3g for LAI, and MTE and FLUXCOM for GPP. As, other than NBP, these reference sets are the same as the ones considered in this paper, the results are also the same. For NBP, despite the different data set, the performance of the models is very similar to the one found for CAMS, our alternative data set. The ILAMB benchmark used in Chapter 5 of the IPCC AR6 (Canadell et al., 2021) also comes to the conclusion of the model having improved from CMIP5 to CMIP6. No qualitative difference can be found between models that have both emission-driven and concentration-driven simulations compared to models with only concentration-driven simulations, and models with both simulations have similar RMSDs in both. This indicates that carbon exchanges are well simulated in these models as the freely evolving fluxes are comparable to results with prescribed atmospheric concentrations.

Centered pattern correlations for these variables and regions are shown in Fig. 17, with a score of 1 meaning the perfect similarity of a model to the reference data, while a value of 0 signifies no relationship. For GPP and LAI, the reference data sets show a very good similarity of above 0.9, while for NBP the differences in the references highlighted in Fig. 12 are highlighted through a small correlation of up to 0.3, with a high anti-correlation in the Southern Hemisphere. Due to this, the precise value of the correlation coefficient between models and reference data set is not a good measure, but it can be seen that the models show a large spread. For GPP, the CMIP6 performance in the tropics is vastly improved, with even higher correlation values for Ncycle models. Other than in the NH, the CMIP5 models show a large spread in correlation values, which has reduced for CMIP6. The correlation distribution for LAI is similar to GPP, with the highest correlation values found in the tropics and globally, but the difference between Ncycle and non-Ncycle models is not as prominent. These overall performance plots underline the specific conclusions from the separate sections above.

Figure 17 Centered pattern correlations between models and reference data sets for annual mean climatology for concentration-driven (a) and emission-driven (b) CMIP5 and CMIP6 models split into Ncycle and non-Ncycle models. A score of 1 denotes the perfect similarity of a model to the reference data, while a value of 0 signifies no relationship. The longer lines denote the MMM, while the grey circle shows the similarity of the alternate data set to the main reference data set. Main and alternate observations are the same as in Fig. 15.

4 Summary and conclusion

To be able to have confidence in future model projections of climate change, Earth system models first need to show the ability to simulate observed climatologies and trends of the carbon cycle in the present-day climate. In the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012), several weaknesses of the simulated carbon cycle were found, such as a general overestimation of photosynthesis and a wide range of values for carbon stocks, which became one of the main areas of focus for improvement for some model groups (Delire et al., 2020). In this study, we have analyzed the land carbon cycle of models participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016a) to investigate whether these weaknesses were improved in the newer model generation, with a special focus on differences arising due to inclusion of an interactive terrestrial nitrogen cycle in some of the CMIP6 models. Concentration- and emission-driven simulations from CMIP5 and CMIP6 models were compared to reference data sets, with 2 out of 18 CMIP5 models and 15 out of 23 CMIP6 models including carbon–nitrogen interactions. We assessed the means, trends, and seasonal cycles of the leaf area index (LAI), the gross primary productivity (GPP), and the land–atmosphere carbon flux (NBP). We furthermore looked at land carbon stocks to see if the large range of values simulated in CMIP5 was reduced in CMIP6.

In general, CMIP6 models show a better performance across all assessed land carbon cycle variables to differing degrees, and no significant differences between the concentration-driven and emission-driven simulation were found in the considered variables that cannot be explained by the number of different models. While there is a bias towards the CLM land component in the CMIP6 models, the different versions (4, 4.5, and 5) do not perform the same, and thus these versions can be seen as independent components for the multi-model mean.

The leaf area index was a weakness of the CMIP5 simulation, as its seasonal cycle was not well captured, and its absolute value was generally overestimated. While the peak of the climatological seasonal cycle of LAI is much better reproduced in CMIP6, the amplitude of the seasonal cycle is weaker in CMIP6 compared to observations due to a weaker drawdown in winter. Thus, LAI should remain an area of focus for future model development. Mean LAI is much better reproduced in CMIP6, while the range of trends in the observations is large enough to cover most models for both CMIP5 and CMIP6. It should be noted that due to correlations between parameters, there are often tradeoffs for the better reproduction of one variable. In CLM5, such a tradeoff had to be weighed between biases for GPP and LAI against high plant functional type (PFT) survivability rates (Lawrence et al., 2019). Therefore, looking at one variable separately instead of the whole model performance can lead to wrong conclusions about the model's ability of reproducing the carbon cycle, depending on which choices were made in the tuning. Similarly, models interactively simulating vegetation cover perform worse in the evaluation of present-day LAI compared to models using observationally derived land cover maps due to simulating trees and grasses in the wrong areas. However, paleo-observations show that large changes in climate, such as during the Paleocene–Eocene Thermal Maximum, have a significant impact on vegetation and its global distribution (e.g., McInerney and Wing, 2011). Therefore, only models with dynamic instead of prescribed vegetation can self-consistently account for future changes in vegetation and the impact of these changes on climate and carbon processes in future projections.

One of the largest improvements due to the inclusion of an interactive nitrogen cycle was seen in GPP, where the CMIP6 nitrogen cycle models were able to capture the seasonal cycle in the Northern Hemisphere well, which was previously overestimated. Beside the improvements in the NH, bias patterns in the tropics showing larger GPP overestimations bordering tropical rainforests are reduced in CMIP6 models, with some of these biases wholly removed in the multi-model mean of the CMIP6 models with interactive nitrogen cycle.

The land carbon sink is underestimated in the Northern Hemisphere, regardless of the CMIP phase or the inclusion of the nitrogen cycle. The models compensate for this by simulating a larger carbon sink in the tropics and the Southern Hemisphere for a global average close to the observed value. An improvement is seen in CMIP6 in capturing the amplitude of the seasonal cycle, which is controlled by carbon uptake through photosynthesis in the growth season and carbon release by respiration. This improvement can largely be attributed to the improved seasonal cycle of GPP.

The large range of soil and vegetation carbon was another large weakness of CMIP5, with inter-model differences of 900 PgC for vegetation carbon and 2500 PgC for soil carbon. This range has not significantly decreased in CMIP6, and it remains an area for improvement.

While we find a significant improvement through the inclusion of the nitrogen cycle for photosynthesis, the effects are reduced for the leaf area index and the land–atmosphere carbon flux. Despite similar NBP values for models with and without interactive nitrogen cycle, models without interactive nitrogen overestimate carbon fertilization, leading to large differences in atmospheric carbon content for future scenario simulations (Kou-Giesbrecht and Arora, 2023). Therefore, the inclusion of further limiting nutrients like phosphorus is important, as they will likely have substantial impacts on future carbon uptake (Yang et al., 2023). Model performance overall has improved from CMIP5 to CMIP6, even with the added complexity introducing more degrees of freedom into the models, as also found in the latest IPCC report (Eyring et al., 2021; Canadell et al., 2021). This is a positive outlook for the future, as many aspects have to be considered when increasing model complexity, such as a need to adjust the existing parameterizations after model structural changes from carbon-only to carbon–nitrogen coupling. Without such adjustments, lowered ecosystem carbon storage simulated by models with N processes would lead to an underestimation of carbon pools (Wang et al., 2022). The increased overall model performance confirms results from the individual model groups that found improved performance in carbon cycle variables compared to previous model configurations, with the biggest improvements seen in LAI and GPP (Ziehn et al., 2017; Danabasoglu et al., 2020). Many areas requiring improvement remain, such as simulated carbon stocks, which saw no significant reduction in the simulated range between CMIP5 and CMIP6, or the inclusion of more nutrient limitations like an interactive phosphorus cycle. The improvement of the carbon cycle in the models since CMIP5 is a step in the right direction for a better understanding and a more accurate simulation of future trends. Based on our analysis, due to the small differences between historical concentration- and emission-driven simulations, despite the increased process realism, we recommend ESMs in future CMIP phases to be based on emission-driven simulations to fully account for climate–carbon feedbacks in future projections, supporting the message from Sanderson et al. (2024). Similarly, due to the significant improvements in GPP with the inclusion of an interactive nitrogen cycle and no detrimental change in the present-day evaluation of any carbon cycle variable, we suggest that the nitrogen cycle should be seen as a necessary part of carbon cycle models in the future.

Appendix A: Carbon cycle in the CMIP models

While a summary of the model components for each model can be found in Tables 1 and 2, some special characteristics of the carbon cycle in the land surface components are given below. As several ESMs use the same land model, the sections are listed by the land model instead of the ESM. If different versions of the same component were used in different CMIP models, then the difference between the versions, as well as the models using each version, are explained.

A1 CABLE + CASA-CNP

The Community Atmosphere–Biosphere Land Exchange model (CABLE; Kowalczyk et al., 2013) version 2.4 is a land surface model coupled to the biogeochemistry module Carnegie–Ames–Stanford Approach carbon cycle model with nitrogen and phosphorus cycles (CASA-CNP; Wang et al., 2010) used in ACCESS-ESM1-5 (Ziehn et al., 2020). CASA-CNP and thus ACCESS-ESM1.5 is the only model in this study to include a phosphorus cycle coupled to the land carbon–nitrogen cycle. A sensitivity study of allowable emissions to nutrient limitation found a reduction in the land carbon uptake by 35 %–40 % with nitrogen limitation and a further 20 %–30 % reduction with nitrogen and phosphorus limitation on the carbon cycle (Zhang et al., 2013), showing the importance of nutrient limitation.

A simple land use scheme accounts for the annual net change in the vegetation tile fractions of each grid cell, which consider 10 vegetated and 3 non-vegetated surfaces. Three live and six dead carbon pools are modeled. Leaf area index (LAI) in ACCESS-ESM1-5 is calculated from a specific leaf area and the size of the leaf carbon pool, while phenology is prescribed. In the CMIP5 model ACCESS-ESM1, which is not considered in this paper due to a lack of variables on the ESGF, LAI was significantly higher than observations, mainly due to an overestimation of LAI in the Northern Hemisphere, despite a significant underestimation of LAI in the tropics. To better match the observations, two parameters were adjusted for ACCESS-ESM1-5, namely one PFT-specific parameter used in the parameterization for the maximum carboxylation rate, and thus related to the nitrogen cycle, and one global parameter related to the daytime leaf respiration rate. Further changes to the model since CMIP5 include the conservation of land carbon, which was not conserved in CMIP5, and the inclusion of wetland tiles in the biogeochemistry calculation and the removal of a spin-up condition which ensured a minimum nitrogen and phosphorus level in soil pools.

A2 CLASS + CTEM

The land component in the Canadian Earth System Model (CanESM) is divided into the physical part represented by the Canadian Land Surface Scheme (CLASS; Verseghy, 1991, 2000; Verseghy et al., 1993) and the biogeochemical processes as simulated by the Canadian Terrestrial Ecosystem Model (CTEM; Arora, 2003; Arora and Boer, 2003, 2005). In the CMIP5 model (CanESM2; Arora et al., 2011), version 2.7 of CLASS was used, while the CMIP6 models CanESM and CanESM-CanOE (Swart et al., 2019a) employ CLASS v3.6. While neither version includes a nitrogen cycle, a parameter representing terrestrial photosynthesis downregulation is included to simulate the effect of nutrient constraints. This parameter is increased in CanESM5 compared to the previous version, CanESM2, resulting in a higher land carbon uptake in CanESM5. Four PFTs are considered in CLASS, while CTEM increases the number to nine PFTs so that phenology can be simulated prognostically. LAI is dynamically simulated, and three live and two dead carbon pools are considered. Added features since CanESM2 include dynamic wetlands and their diagnostic methane emissions.

A3 CLM

The Community Land Model (CLM; UCAR, 2020) is the most commonly used land model in this study, with 11 models across CMIP5 and CMIP6 using 3 different versions of it. For the CMIP5 models, CLM3.5 (Oleson et al., 2008) was used by FIO-ESM, and CLM4 (Lawrence et al., 2011) was used in CESM1-BGC and NorESM1-ME. In CMIP6, CLM4 is used for SAM0-UNICON and TaiESM, which mainly adapted the CESM1 configuration (Lee et al., 2020). CMCC-CM2-SR5 and CMCC-ESM2 use CLM4.5 (Koven et al., 2013), while the newest version, CLM5 (Lawrence et al., 2019), is part of CESM2, CESM2-WACCM, NorESM2-LM, and NorESM2-MM.

While in CLM3.5 the nitrogen limitation was merely represented by a downregulation factor, CLM4 introduced the coupled carbon–nitrogen cycle. Further improvements in CLM4 included transient land cover change modeling, changes to the PFT distribution, and more realistic modeling of permafrost regions. To reduce biases found in CLM4 such as low soil carbon stocks and unrealistic values for GPP and LAI in several regions such as a stark overestimation in the tropics, several parameterizations were changed in CLM4.5. Modifications were made to the canopy processes, including co-limitations on photosynthesis and photosynthetic parameters. Newly introduced features included a vertically resolved soil biogeochemistry with vertical mixing of soil carbon and nitrogen and a more realistic distribution of biological fixation over the year. The structure of the litter and soil carbon and nitrogen pools was adapted to the century model, and ¹³C and ¹⁴C carbon isotopes were introduced.

Finally, in CLM5, many major components of the land model were updated, with a focus on a better representation of land use and land cover change, as well as a more mechanistic treatment of key processes. Changes included a stronger soil moisture control on decomposition, the use of ¹³C and ¹⁴C isotopes for crops, and several changes to the nitrogen cycle and its impact on photosynthesis. Flexible plant C : N ratios were introduced to eliminate instantaneous downregulation of photosynthesis, leaf nitrogen was optimized in the form of the leaf use of nitrogen for assimilation (LUNA; Ali et al., 2016) model, and a model handling the fixation and uptake of nitrogen (FUN; Shi et al., 2016) was included. With respect to the land use and land cover aspect of the model, land unit weights are no longer fixed during the simulation, and the transient PFT distribution was updated. CLM5 considers 22 live and 7 dead carbon pools, as well as 22 PFTs.

A4 CoLM + BNU-DGVM

The Common Land Model (CoLM; Dai et al., 2003), which shares an initial version with CLM but was then developed separately, is the land model component for the CMIP5 model BNU-ESM in the CoLM2005 version. CoLM includes a photosynthesis–stomatal conductance model for sunlit and shaded leaves separately. While carbon–nitrogen cycle interactions were included in the model, they were turned off for the CMIP5 simulations due to not being fully evaluated at the time (Ji et al., 2014).

A5 HAL

The land model for the MRI models in both its CMIP5 version MRI-ESM1 and the CMIP6 version MRI-ESM2-0 is the Hydrology, Atmosphere, and Land surface model (HAL; Hosaka, 2011). It consists of three submodels called SiByl (vegetation), with grass and canopy vegetation layers; SNOWA (snow); and SOILA (soil), with 14 soil layers in the CMIP5 experiments.

A6 ISBA-CTRIP

The land component for the CNRM-ESM2-1 model is presented by the Interaction Soil–Biosphere–Atmosphere (ISBA) land surface model and the total runoff integrating pathways (CTRIP) river routing model (Decharme et al., 2019; Delire et al., 2020). ISBA-CTRIP simulates plant physiology, leaf phenology, carbon allocation and turnover, wild fires, and carbon cycling through litter and soil (Séférian et al., 2019). Land use processes are prescribed instead of simulated, while land cover changes are used to represent anthropogenic disturbances. While the model does not include an interactive nitrogen cycle, its effects are included through an artificial downregulation of photosynthesis and a reduced specific leaf area with increasing CO₂ concentration. There are 6 live and 7 dead carbon pools considered, along with 16 PFTs (Gibelin et al., 2008). Changes since the previous version used in CNRM-ESM-1 include improvements to the photosynthetic and autotrophic respiration schemes.

CTRIP includes carbon leaching through the soil and subsequent transport of dissolved organic carbon to the ocean. As chemical species such as dissolved inorganic carbon are not included, the air–water carbon exchange in the river routing model CTRIP cannot be computed. This leads to a carbon cycle which is not fully bounded.

A7 JSBACH

JSBACH is the land component of the MPI-ESM model, with version 3.2 used for MPI-ESM1.2 (Mauritsen et al., 2019). In the previous version, parameters in the model for the decomposition of dead organic matter were tuned to reproduce the historical atmospheric CO₂ concentrations, with soil and litter carbon stocks merely being the result of this tuning. In version 3.2, decomposition is handled by the YASSO model (Tuomi et al., 2011), based on litter and soil data and resulting in no unconstrained parameters. YASSO simulates four fast soil carbon pools and one slow pool. A total of 18 dead carbon pools is considered due to a different application based on the woody and non-woody origins, as well as above- and belowground decomposition. Additionally, 3 live carbon pools (natural vegetation, crops, and pasture) and 13 PFTs are simulated by JSBACH, while permafrost carbon is not considered. The dynamical vegetation component interacts with the land use changes, modifying the land use data set to conform to the JSBACH setup. JSBACH3.2 includes an interactive terrestrial nitrogen cycle (Goll et al., 2017) driven by the nitrogen demand of the carbon cycle. Further adjustments in v3.2 include the change in carbon timescales in wood pools to be PFT-specific.

A8 JULES

The Joint UK Land Environment Simulator (JULES; Best et al., 2011; Clark et al., 2011) is the land model for the CMIP5 models HadGEM2-CC and HadGEM2-ES with the terrestrial carbon cycle, following the Top-down Representation of Interaction of Foliage and Flora Including Dynamics (TRIFFID; Cox, 2001) dynamic vegetation scheme. The CMIP6 model UKESM1-0-LL (Sellar et al., 2019) employs JULES version 5.0.

Improvements in the version used for UKESM1-0-LL include the introduction of nitrogen cycling, as well as developments to plant physiology and functional types and land use. In this model, nitrogen controls biomass and the leaf area index within TRIFFID, thus only indirectly affecting photosynthetic capacity, as well as limiting the decomposition of litter into soil carbon. For better agreement with observations, global total GPP was tuned down through a reduction in the quantum efficiency of photosynthesis. Furthermore, crop and pasture areas were separated, and a harvest carbon flux was introduced. UKESM1 has four soil carbon pools, nine natural PFTs – increased from five in prior versions – and four PFTs for crop and pasture.

A9 LM

The GFDL Land Model (LM; Anderson et al., 2004) is used in the GFDL Earth System Model with the CMIP5 models GFDL-ESM2G and GFDL-ESM2M (Dunne et al., 2012, 2013), using version 2.0, and the CMIP6 model GFDL-ESM4, employing version 4.1 (Dunne et al., 2020). Neither version includes an interactive nitrogen cycle.

Improvements since CMIP5 include updated soil types in the CORPSE model (Sulman et al., 2014, 2019), hydrology, and radiation, as well as the inclusion of a new fire model FINAL (Rabin et al., 2018), with daily computations instead of previously annual figures and a new model for vegetation dynamics through the perfect plasticity approximation (PPA; Weng et al., 2015). LM4.1 includes 6 live carbon pools for leaves, fine roots, heartwood, sapwood, seeds, and nonstructural carbon, and there are 20 vertical soil levels split into separate fast and slow pools and pools for soil microbes and microbial products. Six PFTs are included, representing C₃ grass, C₄ grass, tropical trees, temperate deciduous trees, and cold evergreen trees. Land use is accounted for through annual wood harvesting, crop planting and harvesting, pasture grazing, and newly included rangelands.

A10 LPJ-GUESS

The Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS; Smith et al., 2014) in combination with the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (H-TESSEL; Balsamo et al., 2009) is the land model used in the EC-Earth models EC-Earth3-CC and EC-Earth3-Veg (Döscher et al., 2022). The model described as “CC” additionally includes ocean biogeochemistry (PISCES) and atmospheric composition for CO₂ (TM5), letting it perform CO₂ emission-driven simulations.

H-TESSEL solves the energy and water balance at the land surface, while vegetation types and vegetation coverage are interactively provided by the coupled LPJ-GUESS, which includes an interactive nitrogen cycle. Compared to the common area-based vegetation schemes, the interactive coupling of LPJ-GUESS to an atmospheric model should improve realism on longer timescales (Döscher et al., 2022). LPJ-GUESS includes 10 L pools, seven vegetation carbon pools, and five soil carbon pools. Wildfires, disturbances, and land use change are simulated on a yearly time step and distributed evenly throughout the year to conserve carbon mass. Land use change dynamics are considered together with a crop module (Lindeskog et al., 2013), including five crop functional types. Three types of plant phenology – evergreen, seasonal deciduous, and stress deciduous – are considered, with only the latter two being simulated with an explicit phenological cycle. Seasonal deciduous PFTs have a fixed growing season length of 210 d, while the growing season for stress deciduous PFTs is determined by a threshold for water stress.

A11 MATSIRO + SEIB-DGVM/VISIT-e

The Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO; Takata et al., 2003) is the physical land model for the MIROC-ESM family – MIROC-ESM and MIROC-ESM-CHEM (with coupled atmospheric chemistry) for CMIP5 (Watanabe et al., 2011) and MIROC-ES2L for CMIP6 (Hajima et al., 2020b) – which consists of a single-canopy layer, three snow layers, and six soil layers down to a depth of 14 m. For the CMIP6 version, a physically based parameterization for snow distribution and snow-derived wetlands was added. Biogeochemistry in MIROC-ESM and MIROC-ESM-CHEM is simulated by the spatially explicit individual-based dynamic global vegetation model (SEIB-DGVM; Sato et al., 2007). It includes 13 PFTs split into 2 for grass and 11 for trees, as well as 2 organic carbon pools. Light capture competition among trees is explicitly modeled instead of parameterized.

MIROC-ES2L for CMIP6 uses the Vegetation Integrative SImulator for Trace gases model (VISIT; Ito and Inatomi, 2012), with changes for coupling to the ESM (adding the “-e” suffix), such as including leaf nitrogen concentrations, and thus limitations to enable fully coupled climate–carbon–nitrogen projections, and land use change processes to get more use out of new land use change (LUC) forcing data sets, such as using five types of land cover. The model does not simulate explicit dynamic vegetation. Three vegetation carbon pools (leaf, stem, and root) are dynamically regulated and have constant turnover rates to three litter and three soil pools. In total, 12 vegetation types are considered. A daily time step is used for the land ecosystem and land biogeochemistry.

A12 ORCHIDEE

The Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE; Krinner et al., 2005; Boucher et al., 2020) land model is used in the IPSL models, with version 1 for the CMIP5 models IPSL-CM5A-LR and IPSL-CM5A-MR and version 2 in the CMIP6 model IPSL-CM5B-LR. The model considers 15 PFTs, as well as 8 vegetation carbon, 4 L carbon, and 3 soil carbon pools. Plant and soil carbon fluxes are computed every 15 min (the same as the atmospheric physics time step), while slow processes like soil and litter carbon dynamics are computed daily instead. The CMIP5 model used a 2-layer bucket model for its soil hydrology, while in CMIP6 an 11-layer soil hydrology scheme is employed. Photosynthesis is parameterized based on the common Farquhar and Collatz schemes for C₃ and C₄, respectively. Nutrient limitation in CMIP6 is introduced through downregulation using a logarithmic function of the CO₂ concentration.

A13 Other: INMCM

The carbon cycle module for INMCM (Volodin, 2007) includes a single soil carbon pool. The most important changes with respect to inmcm4 for the CMIP6 models lie in the atmospheric component of the model, as well as some upgrades to the oceanic component, but there are no changes to the carbon cycle (Volodin et al., 2017b).