Investigating carbon and nitrogen conservation in reported CMIP6 Earth system model data

Reliable, robust, and consistent data are essential foundations for analysis of carbon cycle feedbacks. Here, we consider the data from multiple Earth system models (ESMs) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6). We identify a mass conservation issue in the reported carbon and nitrogen data, with a few exceptions for specific models and reporting levels. The accumulated mass imbalance in the reported data can amount to hundreds of gigatons of carbon or nitrogen by the end of the simulated period, largely exceeding the total carbon–nitrogen pool size changes over the same period. Nitrogen mass imbalance is evident across all reported organic and inorganic pools, with mineral nitrogen exhibiting the most significant cumulative mass imbalance. Due to a lack of detail in the reported data, we cannot uniquely identify the cause of this imbalance. However, we postulate that the carbon mass imbalance primarily arises from missing fluxes in the reported data and inconsistencies between these data and the definitions provided by the C4MIP protocol (e.g., land-use and fire emissions), rather than from an underlying mass conservation issue in the models themselves. Our findings suggest that future CMIP reporting protocols should consider incorporating mass conservation into their data validation processes so that such issues are caught before users have to deal with them, rather than forcing all users to handle this issue in their own way. In addition, attention from model groups to the detailed diagnostic request and definitions, along with their own quality control, will also help to avoid such issues in future. Given that no additional CMIP6 data are currently being published and none are expected in the future, we recommend that data users that rely on a closed carbon–nitrogen cycle address potential flux imbalances by using the workarounds provided in this study.

1 Introduction

The Coupled Model Intercomparison Project Phase 6 (CMIP6, the latest phase at the time of writing) is a collaborative effort among international climate research communities and plays a fundamental role in advancing our understanding of past, present, and future climate (Meehl et al., 2020; Tokarska et al., 2020; Zelinka et al., 2020). CMIP protocols coordinate and standardize climate model experiments, establish consistent requirements for model outputs, and define protocols for model evaluations (Eyring et al., 2016). The Earth system model (ESM) outputs serve as a primary model data resource for climate scientists worldwide, supporting a wide range of analyses (Tang et al., 2024b; Nicholls et al., 2022; Turnock et al., 2020; Stouffer et al., 2017) and forming the scientific basis for assessments such as those conducted by the Intergovernmental Panel on Climate Change (IPCC) to inform global climate policy and decision-making (IPCC, 2021b; Arias et al., 2021). Ensuring the reliability, reproducibility, and robustness of the published data is of paramount importance for advancing our understanding of the climate system.

The Coupled Climate–Carbon Cycle Model Intercomparison Project (C4MIP) is one of the CMIP6-Endorsed MIPs focusing on the design, documentation, and analysis of carbon cycle feedbacks and interactions in climate simulations (Friedlingstein et al., 2006; Jones et al., 2016). The interactions between the carbon cycle and atmospheric CO₂ concentration, as well as physical climate, are some of the most important feedbacks in the earth system (Friedlingstein et al., 2006; Arora et al., 2013, 2020). The nitrogen cycle plays a crucial role in carbon cycle feedbacks. For example, the carbon–climate feedback can become negative when the nitrogen cycle is incorporated into models (Zaehle and Dalmonech, 2011; Zaehle et al., 2015). Observational and modeling studies suggest that nitrogen limitation reduces net primary productivity, a key carbon flux that removes CO₂ from the atmosphere (LeBauer and Treseder, 2008; Schulte-Uebbing and De Vries, 2018). Accurate and robust estimates of carbon cycle feedbacks rely on careful analysis of model output. The model output is one line of evidence but is generally of high importance for the question of carbon cycle feedbacks given that there is, to date, little observational evidence with which to constrain these nonlinear feedback processes. Accurate and robust estimates of carbon cycle feedbacks are also highly desirable for developing climate mitigation and adaptation strategies. The need for robust analysis is why C4MIP made considerable effort to define precise reporting instructions and variable definitions in a way that can be applied to all models (Fig. 1). Tier-1 variables are expected from all models, as is carbon conservation in the reported output (Jones et al., 2016). However, the conservation of carbon and nitrogen in the reported output has not, to date, been assessed.

Figure 1 Land carbon–nitrogen pools and fluxes requested by the C4MIP protocol. The boxes represent pools, while the arrows indicate fluxes and their directions. Details of variable descriptions are provided in Table A1. This figure is reproduced from the C4MIP protocol (Jones et al., 2016) under the Creative Commons Attribution 3.0 License.

The “data request” in CMIP defines specific requirements for model outputs (variables) when submitting experiment results (Juckes et al., 2020). As one of the most complex elements in the infrastructure, the CMIP6 data request encompasses thousands of variables from the various CMIP experiments (Eyring et al., 2016). Despite detailed standardizations and requirements in variable definitions and output formats, the large number of variables presents challenges for data quality control for both modeling groups and the Earth System Grid Federation (ESGF) (Petrie et al., 2021). The data quality control issues lead to challenges for data processing and analyses for research communities who use the data, such as end users who require adjustment of spatiotemporal resolution for various analytical purposes (e.g., estimating carbon cycle feedbacks). Where issues are present, users will either use data that have problems they are unaware of or, rather than being able to use the data directly, have to apply a number of fixes before performing their own analysis.

In this study, we focus on the reported land carbon and nitrogen pools and fluxes. Sections 2 and 3 outline our data collection and processing approach, as well as the method for validating mass conservation. Sections 4 and 5 address mass conservation for land carbon and nitrogen, respectively, including their subpools. Each section begins with an introduction to the theoretical mass conservation equations based on the C4MIP protocol (Fig. 1), followed by an analysis of diagnosed mass conservation issues in the reported CMIP6 carbon–nitrogen cycle data. Section 6 discusses additional issues arising from the CMIP6 data request and C4MIP, which contribute to potential confusions. In Sect. 7, we examine possible causes for mass conservation discrepancies and their implications. Based on the results, we provide suggestions in Sect. 8 for the CMIP data request team, ESM groups, the C4MIP and related MIPs, and CMIP6 data users. We hope this paper will help shed light on these issues, aiding the community in achieving a unified understanding of definitions; reporting protocols and diagnostic recipes; and ultimately providing a cleaner, easier-to-use set of data in CMIP7 to ensure that the remarkable progress made by modeling teams in recent decades is not lost in data reporting issues (Arora et al., 2020, 2013; Jones et al., 2016; Friedlingstein et al., 2006).

2 Data collection and processing

For the purpose of mass conservation validation, the selection of ESMs and experiments is based on the following criteria:

1. The model must provide a relatively complete set of carbon cycle variables requested by C4MIP. This completeness primarily considers the availability of carbon pool sizes; key fluxes such as net primary production (npp) and heterotrophic respiration (rh); and cross-verification fluxes, for instance, net biospheric production (nbp) and netAtmosToLandCO2Flux.

2. The model must have completed three representative experiments: 1pctCO2, historical, and ssp585.

Based on these criteria, this study includes 15 sets of model runs from 12 ESMs across nine institutes. Details on the models and experiments are provided in Table A1.

The C4MIP variables and related data (Fig. 1 and Table A1) were obtained from the ESGF (https://esgf.llnl.gov/, last access: 25 June 2023; data availability varies across different ESMs and experiments, Fig. A1). To ensure consistent data post-processing, the most frequently used monthly data with the model's native grid (with “gn” as “grid_label”; see CMIP6 global attributes details at https://goo.gl/v1drZl, last access: 25 June 2023) were collected.

The data were processed into global-mean, annual-mean values. As a first step, we used the model-specific grid area (“areacella”) and land area fraction (“sftlf”) to determine global land pool sizes and fluxes. To calculate their annual means, we accounted for differences in model-specific calendars by weighting the monthly global pool sizes/fluxes according to the number of days in each month. While the translation from monthly to annual values is not strictly necessary for the mass conservation test, it largely aids in clarity and visualization.

3 Methods

In this paper, mass conservation is defined as follows: for any single pool – whether the total land carbon–nitrogen or their respective subpools – mass is conserved when the net flux (the difference between fluxes entering and leaving the pool) equals the change in pool size over time (Fig. 1). This definition applies at all time points, meaning conservation must hold at each individual time step as well as over longer periods.

To verify mass conservation, we performed the following calculations: first, we reconstructed the net flux by summing the individual influxes and outfluxes for each pool (Fig. 1). Second, we determined the change in the pool size over time. While the first calculation is straightforward, the second one requires consideration. Both the pool size data and flux data are temporal mean values (a requirement from the CMIP6 data request). However, this reporting convention means that there is only one pool size in each month, rather than information about the size of the pool at the beginning and end of the month (which is what you would need to unambiguously calculate the change in the pool size over the month of interest). As a result, blindly calculating the difference in pool size between each reported data point will inevitably result in a net flux that differs from the flux reported in the time series. Considering this, we employed a gradient method (Fornberg, 1988), which calculates gradients using second-order accurate central differences in the interior points and first-order accurate one-sides differences at the boundaries, to derive the change in the pool size from the reported pool size data. It should be noted that this approach is simply a method to obtain a better idea of the actual differential (example provided in Sect. A1 and Fig. A2). It does not fundamentally resolve the discrepancy between differencing discrete data and differentiating continuous data. However, given the relatively long time series of monthly data (e.g., 165 years for the historical period), the numerical errors introduced by the reporting convention and the method used to calculate changes in pool size are expected to be minimal. In other words, while the discrepancy persists, the magnitudes of $\mathrm{d}y/\mathrm{d}t$ derived from both methods should closely align.

4 Carbon mass conservation from reported data

To facilitate comparison and formulation, the reconstructed net flux derived from component fluxes and the reconstructed pool size calculated from subpool sizes are denoted with an asterisk (“*”) for both carbon and nitrogen. The mass imbalance discussed in the following sections, pertaining to both carbon and nitrogen, refers to the discrepancy between the fluxes and the changes in pool sizes within the reported data. The terms “overestimation” and “underestimation” in this paper specifically denote cases where the net flux exceeds or is smaller than the pool size change in the reported data, rather than implying a comparison with observational data or any external benchmark.

4.1 The top-level closure for land carbon: equations and results

Based on the C4MIP protocol and CMIP6 data request (Fig. 1), the Tier-1 land carbon pools and directional composite fluxes should have the following relationships.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\text{cLand}}^{\ast }=\text{cVeg}+\text{cLitter}+\text{cSoil}+\text{cProduct}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{\text{cLand}}^{\ast }}{\mathrm{d}t}}=\text{netAtmosLandCO2Flux}-\text{fCLandToOcean}\end{array}$$

Note the following:

• a.

  We use the sum of the subpools (cLand*) rather than directly using reported cLand in Eq. (2) to ensure consistency in carbon conservation analysis between total land carbon and subpool carbon (cLand and cLand* differ in some ESMs; see details in Sect. 6.1).

• b.

  netAtmosLandCO2Flux and nbp (introduced in CMIP5) are interchangeable based on their definitions (which may also be part of the confusion; see details in Sect. 6.2).

• c.

  fCLeach and fLandToOcean are in the C4MIP protocol but not in the CMIP6 data request and are thus not included in Eq. (2) (see details in Sect. 6.3).

• d.

  fCLandToOcean is only provided by CESM2-WACCM (Fig. A1) with negligible values (on the order of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$, not shown).

Therefore, the following results compare nbp, netAtmosLandCO2Flux, and pool size changes derived from pool size data ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$).

Figure 2 Comparison of net biospheric production (nbp), netAtmosLandCO2Flux, net land carbon flux (fCnetLand*, reconstructed from flux data), and the change of land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

The nbp and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ exhibit nearly identical patterns across all the studied ESMs and experiments (Figs. 2 and A3), with the small differences likely resulting from the differentiation method applied (see details in Sect. 3).

The netAtmosLandCO2Flux and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ either are opposite to each other (the two CMCC models seem to report with a sign issue, Fig. 2) or exhibit the same trend (all the other models with netAtmosLandCO2Flux provided, Fig. 2). According to its definition, the netAtmosLandCO2Flux should be identical to nbp, as demonstrated in all ESMs except for the two CMCC models (Fig. A4). The opposite trends in CMCC models are because of the conflicting description in the C4MIP protocol and the CMIP6 data request (see details in Sect. 6.1). Nonetheless, the consistency among nbp, netAtmosLandCO2Flux (ignoring the sign difference for the CMCC models), and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ indicates three things:

• a.

  Carbon mass should be conserved in the models themselves.

• b.

  The reliability and robustness of the pool size data and aggregated fluxes (nbp and netAtmosLandCO2Flux) are proven.

• c.

  The modeling teams can report data that conserve mass; hence the root of the issue is very likely in confusion about how the reporting protocol is meant to be applied.

4.2 The closure of land carbon from sum of fluxes: equations and results

According to the C4MIP protocol and the CMIP6 data requirements (Fig. 1), the Tier-1 land carbon pools and component fluxes for the net land flux (fCnetLand*) should have the following relationships.

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fCnetLand}}^{\ast }& {\displaystyle }=\text{npp}-\text{rh}-\text{fAnthDisturb}\\ {\displaystyle }& {\displaystyle }-\text{fProductDecomp}-\text{fFireNat}\\ {\displaystyle }& {\displaystyle }-\text{fCLandToOcean}\end{array}\text{(4)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{cLand}}^{\ast }}{\mathrm{d}t}}={\text{fCnetLand}}^{\ast }\end{array}$$

Note that many ESMs do not report fAnthDisturb (Fig. A1). In such cases, the Tier-2 fluxes fHarvestToAtmos and fDeforestToAtmos (Fig. 1) are combined to reconstruct it.

The reconstructed fCnetLand* and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ are well aligned in ACCESS-ESM1-5, CanESM5-1, CanESM5, MIROC-ES2L, and UKESM1-0-LL (Fig. 2), with the difference fluctuating around zero and showing no discernible trends (Fig. A5), suggesting that these differences are due to numerical processing issues alone. In contrast, the reconstructed fCnetLand* is much higher than the change in the pools in all the other studied ESMs (with the green line above the dashed black line, Fig. 2). In these models, the disparity between $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ and fCnetLand* ranges from 2 to 6 Gt C yr⁻¹, with an increasing trend in flux differences, particularly in the 1pctCO2 and ssp585 scenarios (Figs. 2 and A5). The biased imbalance flux indicates that some component fluxes (outfluxes) in Eq. (3) may be missing from the reported data or incorrectly reported/double-counted.

4.3 The significance of cumulative carbon imbalance

Comparing nbp with $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ reveals that, although the differences are relatively small and fluctuate around zero (Figs. 2 and A3), the accumulated mass imbalance is still substantial in certain experiments from some ESMs (Fig. 3). For example, in the 1pctCO2 runs of CESM2-WACCM, NorESM2-LM, and NorESM2-MM, there is an overestimation of approximately 60 Gt C. In contrast, during the historical experiment, nbp and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ are well aligned, resulting in a negligible accumulated imbalance. The mass imbalance in the ssp585 run is also small except for the TaiESM1, where the integration of ($\text{nbp}-\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$) reveals an underestimation of 53.8 Gt C of carbon pool size change, comparable in magnitude to the overestimation resulting from the (${\text{fCnetLand}}^{\ast }-\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$) integration (48.3 Gt C). This indicates that both nbp and fCnetLand* in TaiESM1's ssp585 run are inconsistent with its reported pool size data.

Figure 3 The carbon mass imbalance from the mismatch of net land carbon flux (fCnetLand*, reconstructed from flux data) or net biospheric production (nbp) and the change in land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) for CMIP6 ESMs across different scenarios.

The imbalance flux from the reconstructed fCnetLand* results in an accumulated mass imbalance ranging from 93.7 to 531.8 Gt C, 115.1 to 412.9 Gt C, and 48.3 to 351.2 Gt C for the 1pctCO2, historical, and ssp585 runs, respectively, in models where carbon conservation is not strictly maintained. These values are substantially higher than the cumulative imbalance from $\text{nbp}-\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ (with a maximum of approximately 60 Gt C). By contrast, models that maintain relative carbon conservation (ACCESS-ESM1-5, CanESM5-1, CanESM5, MIROC-ES2L, and UKESM1-0-LL) exhibit much smaller imbalances (<5.5 Gt C), likely reflecting unavoidable numerical errors arising from the differentiation of discrete pool size time series.

From the historical period to ssp585, the land carbon storage trends based on pool size data show a decrease from 1850 to 1975, followed by either continued increase or stabilization until 2100, depending on the ESM (Fig. A6). During the historical period, most ESMs exhibit a decrease in land carbon pool size ranging from 6.8 Gt C (CanESM5-1) to 62.6 Gt C (TaiESM1). In contrast, ACCESS-ESM1-5, CMCC-CM2-SR5, CMCC-ESM2, and MIROC-ES2L show accumulation of land carbon ranging from 39.0 to 109.1 Gt C. Integration of fCnetLand* indicates a maximum carbon accumulation of 465.4 Gt C during the historical period, significantly exceeding the changes observed in land carbon pool sizes (Fig. A6). Additionally, for ESMs exhibiting mass imbalance, none capture the decreasing trend in land carbon before 1975 through integration of fCnetLand*. Under the ssp585 scenario, land carbon pool sizes increase from 2.4 Gt C (TaiESM1) to 600.4 Gt C (CanESM5-1). ESMs with mass imbalance show integration of fCnetLand* at least 1.5 times higher than their corresponding changes in pool sizes. The large discrepancy suggests that ESMs with different fCnetLand* and $\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$ may lead to significant issues in carbon budget (here meaning the internal flows within the carbon cycle, not our remaining carbon budget) assessments. Some of these conservation issues may be solved by considering variables outside the C4MIP protocol (see discussions in Sect. 7.2).

4.4 The closure of land carbon in subpools: equations and results

The subpools within the land carbon pool should also maintain mass conservation. Note that the mass conservation of subpools involves some Tier-2 fluxes, such as rhLitter and rhSoil (subfluxes for rh, Fig. 1), which are unreported for many ESMs (Fig. A1). Therefore, we combine the litter and soil pools when writing the mass conservation equations for subpools.

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fCnetVeg}}^{\ast }& {\displaystyle }=\text{npp}-\text{fDeforestToProduct}\\ {\displaystyle }& {\displaystyle }-\text{fHarvestToProduct}-\text{fAnthDisturb}\\ {\displaystyle }& {\displaystyle }-\text{fVegLitter}-\text{fVegSoil}-\text{fVegFire}\end{array}\text{(6)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{cVeg}}^{\ast }}{\mathrm{d}t}}={\text{fCnetVeg}}^{\ast }\text{(7)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fCnetLitterSoil}}^{\ast }& {\displaystyle }=\text{fVegLitter}+\text{fVegSoil}-\text{fLitterFire}\\ {\displaystyle }& {\displaystyle }-\text{rh}-\text{fCLandToOcean}\end{array}\text{(8)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{cLitterSoil}}^{\ast }}{\mathrm{d}t}}={\text{fCnetLitterSoil}}^{\ast }\end{array}$$

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fCnetProduct}}^{\ast }& {\displaystyle }=\text{fDeforestToProduct}\\ {\displaystyle }& {\displaystyle }+\text{fHarvestToProduct}\\ {\displaystyle }& {\displaystyle }-\text{fProductDecomp}\end{array}\text{(10)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{cProduct}}^{\ast }}{\mathrm{d}t}}={\text{fCnetProduct}}^{\ast }\end{array}$$

Note that the data availability of Tier-2 variables hinders the mass conservation analysis of subpools (Fig. A1). We thus slightly adjust the calculation of the above equations (see details in Sect. 6.4).

Figure 4 The carbon mass imbalance from mismatch of net carbon fluxes for vegetation, litter + soil, and product pools (fCnetVeg*, fCnetLitterSoil*, and fCnetProduct*, reconstructed from flux data) and the change in carbon pool size over time for vegetation, litter + soil, and product pools ($\mathrm{d}\text{cVeg}/\mathrm{d}t$, $\mathrm{d}\text{cLitterSoil}/\mathrm{d}t$, and $\mathrm{d}\text{cProduct}/\mathrm{d}t$, derived from pool size data) for CMIP6 ESMs across different scenarios.

When comparing the reconstructed net flux with the pool size change, it is clear that net flux overestimation or underestimation occurs across all subpools for nearly all the models (Fig. 4; details in Sect. A2 and Fig. A7). Notably, all CanESM models and UKESM1-0-LL show conserved carbon in the vegetation pool. The cumulative imbalance, calculated from the integration of imbalance flux, shows that the vegetation carbon pool exhibits the highest imbalance in most cases, with cumulative mass imbalances ranging from ∼200 to 600 Gt C. Exceptions are found in the 1pctCO2 runs of the two NorESM2 models, where the imbalance from the litter + soil pools surpasses that of the vegetation pools. MIROC-ES2L, despite achieving overall mass conservation in the total land carbon pool (Fig. 3), shows an overestimation of the vegetation carbon pool and an underestimation of the product carbon pool across all experiments (Fig. 4). This discrepancy arises from the absence of anthropogenic fluxes fDeforestToProduct and fHarvestToProduct from cVeg to cProduct in MIROC (Figs. 1 and A1). The CanESM models also conserve total land carbon (Fig. 3); however, their net flux integration reveals an underestimation of the litter + soil carbon pool and an overestimation of the product carbon pool in both the historical and ssp585 runs (Fig. 4). The overestimation in the product pool is due to the absence of fProductDecomp and cProduct in these experiments (Fig. A1), with only the influx fDeforestToProduct reported.

5 Nitrogen mass conservation from reported data

There is no reference composite flux like nbp or netAtmosLandCO2Flux for the land nitrogen cycle in the C4MIP protocol or CMIP6 data request (Fig. 1). Thus, we start by directly analyzing the land nitrogen mass conservation based on the sum of fluxes. We also analyze the mass conservation from the organic and inorganic nitrogen pools. However, because of a lack of available data (Fig. A1), the further investigation of mass conservation of the organic nitrogen subpools (nVeg, nLitter, nSoil, and nProduct) is not included in this study. Besides, the CanESM series models do not include a nitrogen cycle and are therefore excluded from the following figures.

5.1 The closure of land nitrogen from sum of fluxes: equations and results

The Tier-1 land nitrogen pools and component fluxes should have the following relationships (based on the C4MIP protocol, Fig. 1).

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{nLand}}^{\ast }& {\displaystyle }=\text{nVeg}+\text{nLitter}+\text{nSoil}+\text{nProduct}\\ {\displaystyle }& {\displaystyle }+\text{nMineral}\end{array}\text{(12)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fNnetLand}}^{\ast }& {\displaystyle }=\text{fBNF}-\text{fNAnthDisturb}+\text{fNdep}\\ {\displaystyle }& {\displaystyle }+\text{fNfert}-\text{fNloss}\end{array}\text{(13)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{nLand}}^{\ast }}{\mathrm{d}t}}={\text{fNnetLand}}^{\ast }\end{array}$$

The organic and inorganic nitrogen pools should maintain their respective mass conservation as follows.

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{\displaystyle }{\text{nOrganic}}^{\ast }={\displaystyle }\text{nVeg}+\text{nLitter}+\text{nSoil}+\text{nProduct}\text{(15)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fNnetOrganic}}^{\ast }& {\displaystyle }=\text{fBNF}-\text{fNAnthDisturb}+\text{fNup}\\ {\displaystyle }& {\displaystyle }-\text{fNnetmin}\end{array}\text{(16)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\text{nOrganic}}^{\ast }}{\mathrm{d}t}}={\text{fNnetOrganic}}^{\ast }\text{(17)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\text{fNnetMineral}}^{\ast }& {\displaystyle }=\text{fNdep}+\text{fNfert}+\text{fNnetmin}-\text{fNup}\\ {\displaystyle }& {\displaystyle }-\text{fNloss}\end{array}\text{(18)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}\text{nMineral}}{\mathrm{d}t}}={\text{fNnetMineral}}^{\ast }\end{array}$$

Note the following:

• a.

  Similar to the carbon cycle data, we use the sum of subpools (nLand*) rather than the reported nLand for mass conservation analysis (nLand and nLand* differ in some ESMs; see details in Sect. 6.1).

• b.

  Because many ESMs do not report fNloss (Fig. A1), we substitute it with Tier-2 variables, gaseous nitrogen loss (fNgas) and nitrogen leaching (fNleach), where necessary.

Figure 5 Comparison of net land nitrogen flux (fNnetLand*, reconstructed from flux data) and the change of land nitrogen pool size over time ($\mathrm{d}{\text{nLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

In principle, nitrogen mass should be conserved in the organic nitrogen pool, the mineral nitrogen pool, and the total land nitrogen pool (Eqs. 11–18). However, there is a significant discrepancy between the reconstructed net flux and pool size change over time across nearly all land nitrogen pools in all ESMs and experiments (Figs. 5 and A8). Notably, only MIROC-ES2L and MPI-ESM1-2-LR show mass conservation from the reconstructed net land nitrogen flux (fNnetLand*) in all three experiments (Fig. 5). The two NorESM models show fNnetLand* matching the land nitrogen pool size change ($\mathrm{d}{\text{nLand}}^{\ast }/\mathrm{d}t$) in their 1pctCO2 runs.

The gap between the reconstructed net flux and pool size change for the land organic and inorganic nitrogen pools is significantly larger than that for the total land nitrogen pool (Fig. A8, note the larger y axis scale). For the organic nitrogen pool, the fNnetOrganic* and $\mathrm{d}{\text{nOrganic}}^{\ast }/\mathrm{d}t$ only match well in the three experiments from UKESM1-0-LL. For the inorganic nitrogen pool, while in some models, like the two CMCC models and UKESM1-0-LL, the fNnetMineral* and dnMineral*/dt time series are very close to each other (Fig. A8), the cumulative imbalance is still significant considering the small mineral nitrogen pool size (Fig. A9).

5.2 The significance of cumulative nitrogen imbalance

Regarding the cumulative nitrogen imbalance, in ACCESS-ESM1-5, the fluxes underestimate the change in the total land nitrogen and mineral nitrogen pools, while the fluxes overestimate the change in the land organic nitrogen pool, with a mass imbalance ranging from 20 to 80 Gt N (Fig. 6). In CESM2-WACCM, the reconstructed fNnetLand* underestimates the change in the land nitrogen pool size ($\mathrm{d}{\text{nLand}}^{\ast }/\mathrm{d}t$) for the 1pctCO2 and historical experiments by 20.9 and 7.3 Gt N, respectively, while overestimating it by 5.5 Gt N in the ssp585 run. The total land nitrogen changes in CESM2-WACCM for these experimental periods are +11.7, +1.7, and +5.7 Gt N (Fig. A9), respectively. This significant imbalance shows that the change in total land nitrogen varies depending on whether it is assessed using pool size or flux data.

Figure 6 The nitrogen mass imbalance from mismatch of net nitrogen fluxes for total land, land organic, and land inorganic pools (fNnetLand*, fNnetOrganic*, and fNnetMineral*, calculated from flux data) and the change of nitrogen pool size over time for total land, land organic, and land inorganic pools ($\mathrm{d}{\text{nLand}}^{\ast }/\mathrm{d}t$, $\mathrm{d}{\text{nOrganic}}^{\ast }/\mathrm{d}t$, and $\mathrm{d}\text{nMineral}/\mathrm{d}t$, derived from pool size data) for CMIP6 ESMs across different scenarios.

Examining the mass conservation of the organic and mineral land nitrogen subpools, the net flux in CESM2-WACCM significantly underestimates the organic nitrogen pool size change by 136.9 to 273.6 Gt N while overestimating the mineral nitrogen pool size change by 142.4 to 252.7 Gt N (Fig. 6). CMCC-CM2-SR5 and CMCC-ESM2 both overestimate the total land nitrogen and organic nitrogen while underestimating the mineral nitrogen, with a mass imbalance of less than 25 Gt N. For CMCC-CM2-SR5, the fNnetLand* consistently exceeds the land nitrogen pool size change, with overestimation ranging from 12.6 to 19.7 Gt N. This imbalance accounts for at least 5 times the total land nitrogen change for the experimental period in this model (<2.2 Gt N). MIROC-ES2L, MPI-ESM1-2-LR, NorESM2-LM, and NorESM2-MM are relatively mass-conserved in their total land nitrogen pool (mass imbalance <5 Gt N), particularly in their 1pctCO2 runs (mass imbalance <0.8 Gt N). However, this conservation is achieved at the expense of simultaneously underestimating the land organic nitrogen pool size while overestimating the mineral nitrogen pool size (for MIROC-ES2L) or vice versa (for the other three). Both the overestimations and underestimations are substantial, ranging from 25 to 75 Gt N. The reconstructed net flux in TaiESM1 shows an overestimation of its land nitrogen pool size in the historical and ssp585 experiments (∼25 Gt N, with no data provided for its 1pctCO2 run), primarily due to the overestimation of the organic nitrogen pool. In UKESM1-0-LL, the reconstructed net flux overestimates the total land nitrogen by 4.1 Gt N in the 1pctCO2 experiment and underestimates it by 6.7 and 7.0 Gt N in the historical and ssp585 experiments, respectively. The mass imbalance for UKESM1-0-LL predominantly arises from its organic nitrogen pool.

Based on reported pool size data, the total land nitrogen changes across all studied ESMs during the 1pctCO2, historical, and ssp585 experimental periods range from −2.0 to +13.1 Gt N, +0.4 to +3.8 Gt N, and −1.7 to +5.5 Gt N, respectively (Fig. A9). Changes in the mineral nitrogen pool are substantially smaller, remaining below 0.3 Gt N for all experiments and ESMs, with the exception of the 1pctCO2 run in NorESM2-LM (Fig. A9). Comparatively, the cumulative nitrogen imbalance from the reconstructed net flux (Fig. 6) is orders of magnitude larger than the total nitrogen pool size changes over the same periods.

6 Other issues with the reported CMIP6 data

6.1 The reported cLand and nLand

In many of the ESMs examined in this study, the reported total land carbon or nitrogen (cLand or nLand) does not match the sum of the respective subpools (i.e., the reconstructed cLand* or nLand*). The discrepancy averages between 50 and 150 Gt C or 25 and 55 Gt N across the three experimental periods. This inconsistency is unexpected, as pool size, being a state variable, should be simple to calculate and report. This hence indicates that modeling teams may face challenges in fully adhering to C4MIP guidelines, an issue that warrants attention as we approach CMIP7. However, without further evidence (higher-resolution or more comprehensive datasets), we are unable to analyze the underlying causes in the CMIP6 data.

6.2 The definition of nbp and netAtmosLandCO2Flux

The netAtmosLandCO2Flux was a new variable introduced in the CMIP6 phase of C4MIP, defined as the “net flux of CO₂ between atmosphere and land (positive into land) as a result of all processes.” Many models do not report this flux (Fig. A1). Instead, all models report net biospheric production (nbp), defined as the “carbon mass flux out of atmosphere due to net biospheric production on land”, a variable introduced in CMIP5. Based on their definitions, netAtmosLandCO2Flux and nbp appear interchangeable.

However, it should be noted that in the C4MIP paper (Jones et al., 2016), netAtmosLandCO2Flux is illustrated as directed from land to atmosphere (Fig. 1) and described as the “total flux of CO₂ from land to atmosphere,” which conflicts with its definition in the CMIP6 data request and may be misleading (e.g., the opposite signs for netAtmosLandCO2Flux and nbp observed in the CMCC models, Fig. 2). Resolving these issues may help reduce confusion and improve the consistency of reported data.

6.3 The changes in fCLeach and fLandToOcean from C4MIP protocol to CMIP6 data request

The original C4MIP design includes a leaching carbon flux (fCLeach) as a Tier-1 variable (Fig. 1). However, this flux is not in the CMIP6 data request (and therefore the ESGF) and is thus excluded from Eq. (2). Additionally, fLandToOcean in the C4MIP protocol is renamed as fCLandToOcean in the CMIP6 data request. These changes introduce potential confusion for data users and suggest that a living variable definition and protocol may be essential for providing groups with clarity, given the slow-moving nature of paper publication (which makes papers a source that quickly becomes outdated). Living data/protocol papers may provide another solution which could be updated quickly enough to avoid becoming quickly out of date.

6.4 The adjustments for conservation calculations due to limited data availability

For carbon conservation in vegetation and litter-soil pools, consideration of the Tier-2 component fluxes (fVegFire and fLitterFire) for fFireNat is unavoidable (Eqs. 5 and 7). However, fFireNat itself is reported only in TaiESM1 and the two NorESM models (Fig. A1), and component fluxes are even rarer (only CESM2-WACCM reports fVegFire, and none of the studied models report fLitterFire, Fig. A1). Consequently, in calculating fCnetVeg* and fCnetLitterSoil, we replace fVegFire in Eq. (5) with fFireNat and omit fLitterFire from Eq. (7). While this adjustment is not ideal, it is the only feasible approach given the limited data availability.

It is common for ESMs to report Tier-2 variables while omitting Tier-1 variables for certain fluxes (e.g., fHarvestToAtmos and fDeforestToAtmos vs. fAnthDisturb; fNgas and fNleach vs. fNloss, Figs. 1 and A1). This does not pose a major issue, as we can reconstruct the composite variables by combining the component variables, as done in the conservation calculations. However, this suggests that the C4MIP or CMIP6 data request could consider upgrading the priority of some component variables and potentially removing certain composite variables.

7 Discussions

7.1 Top-level implications

First, it is important to note that we cannot directly determine whether there is a true “mass conservation issue” in the physics of the ESMs based on the published data. From a modeler's perspective, such an issue seems unlikely. Furthermore, we cannot directly verify discrepancies that may arise from the post-processing of the original ESM output to meet the CMIP6 data requirements. It is worth noting that CMIP6-requested variables, the most common of which is monthly gridded data, do not match the original temporal solving frequency in ESMs (Thornton et al., 2007; Sokolov et al., 2008; Wiltshire et al., 2021; Lawrence et al., 2019). This means that post-processing steps like temporal averaging and spatial summing/regridding are inevitable. Numerical errors or other issues may occur during this data processing. Without thoroughly examining the ESM groups' post-processing steps for CMIP6 data submission, it is impossible to pinpoint the source of the mass imbalance. Therefore, the following discussion of the mass conservation issue refers solely to the notable mismatch between the flux data and pool size data submitted to CMIP6/the ESGF by ESMs – based on the C4MIP protocol.

The inconsistency between changes in pool size (derived from pool size data) and the reconstructed net flux (derived from flux data, as detailed in Sects. 4 and 5) is substantial and is observed across various models and experiments for both carbon and nitrogen, including their subpools (however limited data availability for subpools precludes further discussion here). This discrepancy warrants careful attention in future data requests and presents challenges for current data users.

The carbon mass conservation analysis reveals that the land carbon pool size data align well with both the nbp and netAtmosLandCO2Flux (where available, and ignoring sign issues) in most cases (Fig. 2), which enhances the credibility of the pool size data. The matching of nbp and pool sizes (state variables) also supports the assertion that the ESMs maintain mass conservation in their raw data.

7.2 The inconsistent land carbon pool size change ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$) and net flux (fCnetLand*)

To understand the potential causes of the land carbon conservation issue, we look into the component fluxes in the net flux (Eq. 3), which include npp, rh, fAnthDisturb, fFireNat, fProductDecomp, and fLandToOcean. For clarity, Table 1 provides a summary of the land carbon mass conservation status across the models and experiments analyzed.

Table 1 Land carbon mass conservation based on balance between reported fluxes and changes in pool sizes (Eqs. 3 and 4).

Y: mass is conserved (cumulative imbalance ≦5 Gt C). N: mass is not conserved (cumulative imbalance >5 Gt C). N (−): mass is not conserved, but including this variable decreased mass imbalance (make things better). N (+): mass is not conserved, and including this variable increased mass imbalance (make things worse). N (n.d.): mass is not conserved, no data for this variable. /: mass is conserved with C4MIP variables, so no further processing.

7.2.1 fAnthDisturb vs. fLuc: the land-use emissions reported outside C4MIP

The fAnthDisturb variable, comprising fDeforestToAtmos and fHarvestToAtmos (Fig. 1), is related to land use, land-use change, and forestry (LULUCF) emissions (Houghton et al., 2012; Gasser and Ciais, 2013; Stocker and Joos, 2015; Lawrence et al., 2016; Hurtt et al., 2020). Among the studied ESMs, only UKESM1-0-LL – showing negligible mass imbalance (Figs. 2 and 3) – reports an increasing fHarvestToAtmos from 1 to approximately 5 Gt C yr⁻¹ during the historical to ssp585 periods. The historical fAnthDisturb in UKESM1-0-LL aligns with modeled LULUCF emissions from various dynamic global vegetation models (DGVMs; process-based) and bookkeeping models (e.g., OSCAR, used in the Global Carbon Budget program) (Gasser et al., 2020).

Except for UKESM1-0-LL, CanESM5, CanESM5-1 (both mass-conserved in the land carbon pool), and CESM2-WACCM (not mass-conserved) have reported fDeforestToAtmos (0–1.5 Gt C yr⁻¹). The other ESMs have not provided any data for fAnthDisturb, fDeforestToAtmos, or fHarvestToAtmos (Fig. A1). The definitions of these fluxes are similar to those of gross land-use emissions. Their absence in the historical and SSP585 experiments appears unrealistic and may contribute to the diagnosed carbon mass conservation issue. This is partially corroborated by the mass conservation analysis of the land carbon subpools (Sect. A2 and Fig. A7), where the vegetation carbon pool – which also has fAnthDisturb (or fDeforestToAtmos, and fHarvestToAtmos) as an outflux (like the total land carbon pool, Fig. 1) – accounts for most of the land carbon mass imbalance in several ESMs (Fig. 4).

To further investigate the land-use fluxes, we analyze the fLuc (“net carbon mass flux into the atmosphere due to land-use change”) variable from the CMIP5 data request. With the exception of ACCESS-ESM1-5, CanESM5-1, CanESM5, and MIROC-ES2L (in which land carbon is mass-conserved, Fig. 3), all other ESMs report fLuc (Fig. A1), with values $\sim \mathrm{4}\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ for the two NorESM models and $<\mathrm{1.5}\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ for most of the others. Since fLuc is not included in the C4MIP protocol (Fig. 1), it does not fit into the equations used for mass conservation calculations (see Sect. 3). Given its definition, we assume fLuc to be an additional outflux from the land carbon pool and recalculate the imbalance accordingly.

NorESM2-MM and TaiESM1 are mass-conserved in their historical runs when fLuc is considered, with carbon mass imbalances dropping from 226.5 and 115.1 Gt C to 4.7 and 3.1 Gt C, respectively (Table 1 and Fig. A10). The mass imbalance in the NorESM2-LM's ssp585 run also significantly decreases from 351.2 to 22.0 Gt C (Fig. A10). Other ESMs with fLuc considered show varying degrees of imbalance reduction (Fig. A10), although the imbalances remain substantial (>55 Gt C). The fLuc results further highlight the importance of including land-use emissions for mass conservation. However, the discrepancy between the C4MIP protocol and ESM data reporting obscures this issue and may cause confusion.

7.2.2 fFireNat vs. fFire: the fire emissions reported outside the C4MIP

NorESM2-LM, NorESM2-MM, and TaiESM1 are the three models that report fFireNat (Fig. A1). However, the mass imbalance in these four models remains substantial (Figs. 2 and 3). We notice there is a similar CMIP5-requested variable with fFireNat, fFire (“carbon mass flux into atmosphere due to CO₂ emission from fire excluding land-use change”), which is reported by many ESMs (Fig. A1). Thus, we re-check the mass conservation by replacing fFireNat with fFire, where fFire is available.

It is found the fFire is quite high in the 1pctCO2 experiment for the models reporting this flux (e.g., in the range of 1 to 4 Gt C yr⁻¹ with an increasing trend). Considering fFire in the 1pctCO2 run leads to the mass conservation of MPI-ESM1-2-LR and TaiESM1, with the imbalance of 1.6 and 9.1 Gt C, respectively (Table 1 and Fig. A11). The imbalance also drops from 432.7 to 59.7 Gt C in CESM2-WACCM and from 531.8 to 207.0 Gt C in CMCC-ESM2 (Fig. A11). In the historical and ssp585 experiments, though the consideration of fFire largely decreases the imbalance in CMCC-ESM2, it results in some negative imbalance in CESM2-WACCM and MPI-ESM1-2-LR (i.e., fFire is too large). That indicates there might be double counting when simply replacing fFire with fFireNat, which complicates the analysis.

7.2.3 The remaining component fluxes

Besides fAnthDisturb and fFireNat, there are four other fluxes determining the net land carbon fluxes (Eq. 3): npp, rh, fProductDecomp, and fCLandToOcean. Among them, npp and rh stand out as the major carbon fluxes, extensively studied and analyzed in numerous investigations (Jian et al., 2022; Wei et al., 2022; Zhu et al., 2022; Qiu et al., 2023). Observations contribute to constraining the uncertainty of modeled npp and rh (Haverd et al., 2013; Guenet et al., 2024), thereby enhancing their data credibility.

The flux fProductDecomp is reported by ACCESS-ESM1-5, CESM2-WACCM, CMCC-CM2-SR5, CMCC-ESM2, MIROC-ES2L, MPI-ESM1-2-LR, and UKESM1-0-LL (Fig. A1). Notably, four of these seven ESMs demonstrate land carbon mass conservation (Fig. 3). During the historical to SSP585 period, fProductDecomp increases significantly, reaching up to 5 Gt C yr⁻¹ in CESM2-WACCM and MPI-ESM1-2-LR and as high as 7 Gt C yr⁻¹ in MIROC-ES2L. This substantial flux underscores its importance in closing the carbon cycle. The product pool mass imbalance observed in MIROC-ES2L highlights the importance of carbon fluxes into the product pool (Fig. 4).

Among the studied ESMs, CESM2-WACCM is the only model that provides fCLandToOcean. As most of the CMIP6 ESMs may not simulate the river transport and we have not found other similar CMIP6 variables representing the carbon from land to ocean, we do not discuss this flux here.

7.3 The complexity of carbon mass conservation for subpools

In principle, land carbon subpools should also adhere to mass conservation (Eqs. 5–10). However, the results indicate that most do not, even for those ESMs where total land carbon is mass-conserved (Figs. 4 and A7). Given that fluxes for subpool mass conservation follow specific pathways – such that an outflux for one subpool serves as an influx for the next (e.g., fVegLitter is an outflux for cVeg and an influx for cLitter) – imbalanced subpools present challenges in pinpointing the potentially problematic fluxes. Therefore, it would be unreasonable to attribute the imbalance to any single flux based solely on the reported data. This complexity necessitates a more detailed examination of the underlying processes and how these fluxes are reported to identify sources of mass imbalance in the subpools (see Sect. 8 below).

7.4 Nitrogen mass conservation: an emerging and challenging issue

The conservation of land nitrogen mass presents a more intricate challenge due to the complexities inherent in the nitrogen cycle (e.g., the need to conserve both organic and mineral nitrogen, Fig. 1). The total land nitrogen, as well as the organic and mineral nitrogen pools, should maintain its respective mass conservation (Eqs. 11–18). However, current CMIP6 data only supports mass conservation in the total land nitrogen pool for MIROC-ES2L, MPI-ESM1-2-LR, NorESM2-LM, and NorESM2-MM and in the organic nitrogen pool for UKESM1-0-LL (Fig. 6). The imbalance in the reported mineral nitrogen pool is also significant across all the studied CMIP6 ESMs (Fig. 6).

Unlike the carbon cycle, there is no reference nitrogen flux equivalent to nbp or netAtmosLandCO2Flux in the CMIP6 data request and C4MIP design, which complicates the direct validation of both nitrogen pool size data and flux data (Fig. 1). The current fNloss includes fluxes into both the atmosphere (fNgas) and the ocean (fNleach), making it challenging to define a clear reference flux (Fig. 1). Additionally, fNup and fNnetmin link the organic and mineral pools in two directions (Fig. 1), further complicating the attribution of the imbalances.

It is worth noting that only MPI-ESM1-2-LR and UKESM1-0-LL have reported the fNAnthDisturb flux. While land nitrogen in MPI-ESM1-2-LR is mass-conserved and in UKESM1-0-LL it is not, this suggests that attributing the total land nitrogen imbalance solely to anthropogenic disturbance may not be appropriate. The lack of a straightforward reference flux and the complex interconnections within the nitrogen cycle highlight the need for more detailed scrutiny to identify and address the sources of nitrogen mass imbalance in these models.

7.5 The quest for reporting balanced carbon and nitrogen cycles

The carbon and nitrogen cycle data, together with many climate variables, are frequently used to support analyses with various purposes (IPCC, 2021a; Arora et al., 2020; Thornhill et al., 2021; Canadell et al., 2021; Hermans et al., 2021; Fan et al., 2020). These results are particularly crucial for reduced complexity models (emulators), decision-making, and climate policy (Meinshausen et al., 2022; Koven et al., 2022; Kikstra et al., 2022). Considering the challenges in data storage and distribution, it is reasonable that CMIP requests the “most useful” data rather than universally useful output. There is no “one solution for all” for the CMIP data request, meaning that different research groups require varying levels of post-processing based on their specific purposes. This means that ensuring consistency between the model's original outputs and the processed data is important.

Ensuring mass conservation is crucial for data users, especially those analyzing the comprehensive behaviors of the full carbon–nitrogen cycle, such as reduced complexity modeling groups (Meinshausen et al., 2011a, b; Nicholls et al., 2021, 2020; Tang et al., 2024b). Estimating carbon cycle feedback also relies on published ESM data (Arora et al., 2020; Melnikova et al., 2021). Although the nitrogen cycle is relatively new in ESMs (Davies-Barnard et al., 2022), its significance has been widely discussed and demonstrated in both experimental results and models (Zaehle et al., 2010, 2014; Zaehle and Dalmonech, 2011; Schulte-Uebbing and De Vries, 2018). Addressing mass conservation of carbon and nitrogen data, a fundamental requirement for accurate analysis, remains a significant challenge. Achieving mass conservation in reported data would significantly simplify the life of users, facilitating more science done faster.

8 Suggestions for the path forward

Based on the preceding discussions, the following practical recommendations are provided for the CMIP data request team, modeling groups, C4MIP and related MIPs (for the future development), and CMIP6 data users, focusing on carbon–nitrogen cycle data.

8.1 For the CMIP data request and ESM groups

• a.

  It is suggested that the next generation of the CMIP data request and the ESM groups incorporate mass conservation of reported data as a data validation routine. ESM groups may consider reporting fluxes that are consistent with the pool size changes; i.e., there may need to be a better way to represent the assumptions required to ensure consistency between the sum of the reported fluxes and the change in pool sizes. Since the process from original model outputs to reported data is not fully transparent to data users, any imbalances in the reported data, traceable to the smallest subpools the model can analyze, must be documented by the ESM groups. This documentation should include the destination of the imbalance flux (e.g., atmosphere or ocean) and the underlying causes (model process or technical issues like temporal averaging). The ESGF may need to develop necessary data quality control tools considering mass conservation for data publication. Given the significant effort required to manage and share large data volumes, it is reasonable to prioritize these fundamental checks to guarantee that the data shared are both reliable and valuable to the broader community.

• b.

  Direct communication between the team behind the CMIP data request (based on MIP protocols) and the ESM groups should be strengthened to address data reporting challenges and to inform updates to MIP protocols accordingly. Currently, incomplete data reporting (Fig. A1) or fluxes that do not adhere to the prescribed definitions (as evidenced by fLuc and fFire) indicate a disconnect between ESMs and MIP protocols, as ESMs may find it challenging to fully align with these protocols. Bridging this gap requires a collaborative effort and should be addressed promptly to improve data quality and make the substantial data-sharing effort done by the ESGF worth it.

• c.

  The CMIP data request may consider mandating complete data submission from ESM groups for MIPs (e.g., a core set of required variables like the Tier-1 variables in C4MIP) before publication. In particular, reporting Tier-1 variables should be a prerequisite for reporting Tier-2 variables to prevent confusion among data users. For variables requested by MIPs that are not simulated by ESMs, ESM groups should provide a data availability statement for the variable rather than simply not reporting them. Specifically for C4MIP, the note may indicate how to conserve the carbon and nitrogen despite the absence of this variable. Such practices would enhance mass conservation validation for both data users and ESM groups. Currently, the absence of certain variables (Fig. A1) leads to ambiguity from the perspective of data users: it is unclear whether these variables are not modeled or simply not reported.

8.2 For C4MIP and related MIPs

• a.

  For the carbon cycle, collaborations between C4MIP, the Land Use Model Intercomparison Project (LUMIP) (Lawrence et al., 2016), the Fire Model Intercomparison Project (FireMIP) (Rabin et al., 2017), and possibly others should be considered to harmonize data requests and avoid overlap in variables with similar definitions. Currently, the coexistence of variables such as fLuc and fAnthDisturb, or fFire and fFireNat, not only places an additional burden on ESM teams for data reporting but also creates confusion for data users. From a user perspective, relying on model-specific variables to achieve mass conservation is far from ideal. Addressing this issue should be a priority in future updates to the MIP protocols.

• b.

  It is suggested to promote fDeforestToAtmos and fHarvestToAtmos (Tier-2 fluxes for fAnthDiturb), fVegFire and fLitterFire (Tier-2 fluxes for fFireNat), and rhLitter and rhSoil (Tier-2 fluxes for rh) to Tier-1 variables as they are necessary for analyzing mass conservation in key subpools.

• c.

  In the nitrogen cycle, it is suggested to add a new Tier-1 variable to represent the net nitrogen flux between land and atmosphere (e.g., netAtmosLandNgasFlux, similar to nbp and netAtmosLandCO2Flux in the carbon cycle). If ESMs calculate this flux internally (thus, no data processing for CMIP data submission involved), it will aid ESM groups and data users in identifying whether any imbalance originates from fluxes or pool size discrepancies.

• d.

  It is also suggested to promote fNgas and fNleach to Tier-1 variables and remove fNloss, considering their distinct destinations (atmosphere vs. ocean). This might be necessary for the future analyses of the land nitrogen cycle and its interaction with the atmosphere and ocean.

8.3 For users of existing CMIP6 data

• a.

  It is recommended to prioritize the use of carbon pool sizes (cVeg, cLitter, cSoil, and cProduct), net primary production (npp), heterotrophic respiration (rh), and net biospheric production (nbp). For instance, in carbon cycle feedback calculations, utilizing consistent carbon pool size differencing, as exemplified in CMIP6 carbon cycle feedback analyses (Arora et al., 2020), while avoiding time-integrated fluxes used in CMIP5 and earlier studies (Arora et al., 2013; Friedlingstein et al., 2006), is advisable.

• b.

  When a closed carbon–nitrogen cycle is strictly required, specifying how imbalance fluxes are calculated and managed becomes crucial. Data users may consider two main approaches for carbon mass conservation: (i) directly using npp and rh flux data and manually calculating the imbalance flux, attributing it to disturbances (land-use or fire emissions, which should not be assumed as zero in many models/experiments) or land-to-ocean carbon transport flux and (ii) using npp, rh, fLuc, and fFire (where available) and manually calculating the imbalance flux, attributing it to other disturbances or land-to-ocean transport.

• c.

  For nitrogen mass conservation, prioritizing the use of pool sizes (nVeg, nLitter, nSoil, nMineral, and nProduct), biological nitrogen fixation (fBNF; closely linked to npp), nitrogen deposition (fNdep; verifiable through the input4MIPs dataset), nitrogen uptake (fNup), and net mineralization (fNnetmin; both significant nitrogen fluxes) is recommended. Then, attributing imbalance fluxes in mineral and organic nitrogen pools to mineral nitrogen loss and anthropogenic perturbation, respectively, might be considered with minimal data adjustments.

9 Conclusions

The CMIP6 data archive serves as the primary data repository for various communities concerned with climate change. Building and maintaining such a vast data infrastructure is a collaborative effort, understandably fraught with challenges. In this study, we have identified a common issue of mass conservation in reported CMIP6 data from ESMs and experiments following the C4MIP protocol, resulting in accumulated mass imbalances reaching hundreds of gigatons of carbon and nitrogen. The potential reasons include incomplete data submission, inconsistent definitions of fluxes with the C4MIP protocol, different definitions of reporting, particularly concerning land-use and fire emissions, and numerical errors from spatiotemporal regridding (minor). The presence of mass imbalance limits the usability of the data and may exacerbate uncertainties in results/conclusions derived from the data. Therefore, we recommend that ESM groups and the next generation of CMIP, C4MIP, and related MIPs prioritize mass conservation as a fundamental data validation step and document any identified mass imbalances. We think this is a reasonable minimum step to justify the effort required to process, share and analyze the data included in CMIP. We also suggest that data users that require closed carbon–nitrogen cycles prioritize carbon–nitrogen pool size data and major carbon–nitrogen flux data (ignoring other data which would conflict with the mass balance inferred from these sources), while also specifying how they handle any remaining imbalance in their analyses (for example, attributing these imbalances to anthropogenic disturbances).

Achieving these goals requires collaborative efforts among the CMIP data request team, related MIPs, and modeling groups. Given the substantial efforts already invested in making ESM outputs publicly accessible and analyzable, conducting mass conservation verification is valuable. This step not only facilitates data analysis and reanalysis but also reduces uncertainty in the results and conclusions.

Appendix A

A1 Methods for calculating the change in the size of the pools

The reported pool size data represent the mean value for each year, which provides no information about the size of the pool at the beginning or the end of the year. Discrete differencing of the pool size data using $\left[\text{pool size}\right(t+\mathrm{1})-\text{pool size}(t\left)\right]$ actually calculates the difference between the midpoints of year (t+1) and year (t), rather than the within-year change. In contrast, the reported flux data represent the average flux over the course of year (t). This discrepancy means that the pool size change derived from discrete differencing $\left[\text{pool size}\right(t+\mathrm{1})-\text{pool size}(t\left)\right]$ does not correspond to the net flux derived from the reported flux data, potentially leading to a biased estimation of the net flux.

The following provides an example to illustrate the difference between gradient differencing and discrete differencing. Assume a constant flux of 1 Gt C yr⁻¹ over 10 years, with the pool size starting at zero. The true pool size time series would be [$\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{10}$] (11 data points, unit: Gt C). However, since the reported pool size is the mean value for each year, the pool size data would be [$\mathrm{0.5},\mathrm{1.5},\mathrm{2.5},\mathrm{\dots },\mathrm{9.5}$] (10 data points, unit: Gt C, Fig. A2). Discrete differencing of this pool size data yields a flux of 0.5 Gt C yr⁻¹ for the first year (pool size data at t=1 minus zero start), which is lower than the actual flux of 1 Gt C yr⁻¹. By applying gradient differencing, this issue can be avoided.

The same principle can be applied to more complex examples. The point here is that this is a clear issue, but it is also a small one because we have relatively long time series and relatively well-behaved data (so, while we do not have information about the exact values at the start and of the year, the other information provides enough information to constrain their possible values well).

Figure A1 Overview of data availability for the studied CMIP6 ESMs and experiments. The green dots indicate the data are available from the ESGF.

Table A1 Details of variables, experiments, and models.

^a: The experiment label is called “variant_label” in CMIP6 global attributes (https://goo.gl/v1drZl, last access: 25 June 2023). It is a label constructed from four indices stored as global attributes: “r” is followed by the realization index, “i” by the initialization index, “p” by the physics index, and “f” by the forcing index.

Figure A2 The comparison of differencing methods for pool size data time series.

Figure A3 The difference (imbalance flux) between net biosphere production (nbp) and the change of land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

Figure A4 The difference (imbalance flux) between netAtmosLandCO2Flux and the change of land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

Figure A5 The difference (imbalance flux) between “net land carbon flux” (fCnetLand*, reconstructed from flux data) and the change of land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

Figure A6 The mismatch of changes in land carbon pool size derived from pool size data (cLand*(t) − cLand*(t0)) (solid lines) and those retrieved from the integration of net land carbon flux reconstructed from flux data (fCnetLand*, dashed lines) in CMIP6 ESMs across different scenarios.

A2 The carbon mass imbalance from the land carbon subpools

For the vegetation carbon pool, the reconstructed net flux (fCnetVeg*) consistently overestimates the vegetation carbon pool size change ($\mathrm{d}\text{cVeg}/\mathrm{d}t$) in all experiments for CESM2-WACCM, CMCC-CM2-SR5, CMCC-ESM2, MIROC-ES2L, and TaiESM1 (Fig. A7). The fCnetVeg* overestimates the pool size change in the historical and ssp585 runs of both NorESM models, while underestimating the change in the 1pctCO2 run. All CanESM models and UKESM1-0-LL show matching fCnetVeg* and $\mathrm{d}\text{cVeg}/\mathrm{d}t$.

Figure A7 Comparison of net carbon flux for vegetation, litter + soil, and product pools (fCnetVeg*, fCnetLitterSoil*, and fCnetProduct*, reconstructed from flux data) and the change of carbon pool size over time for vegetation, litter + soil, and product pools ($\mathrm{d}\text{cVeg}/\mathrm{d}t$, $\mathrm{d}\text{cLitterSoil}/\mathrm{d}t$, and $\mathrm{d}\text{cProduct}/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

The mismatch between net flux and pool size change is less pronounced for the litter + soil carbon pool (fCnetLitterSoil* and $\mathrm{d}\text{cLitterSoil}/\mathrm{d}t$) compared to the vegetation pool (Fig. A7). Nonetheless, CESM2-WACCM, CMCC-CM2-SR5, and CMCC-ESM2 show a significant overestimation by the reconstructed net flux (${\text{fCnetLitterSoil}}^{\ast }>\mathrm{d}\text{cLitterSoil}/\mathrm{d}t$) in all experiments. The reconstructed fCnetLitterSoil* also overestimates the litter + soil pool size change in all experiments from NorESM2-LM and NorESM2-MM, though the historical run shows a smaller overestimation compared to their 1pctCO2 and ssp585 runs (Fig. A7). All CanESM models, MIROC-ES2L, TaiESM1, and UKESM1-0-LL, exhibit closer alignment between fCnetLitterSoil* and $\mathrm{d}\text{cLitterSoil}/\mathrm{d}t$, although cumulative imbalances remain substantial in certain cases (Fig. 4).

The product pool presents a smaller pool size change compared to the vegetation and litter + soil pools (Fig. A7, around zero where the product pool size is reported). The comparison of the reconstructed net flux (fCnetProduct*) and the derived pool size change over time ($\mathrm{d}\text{cProduct}/\mathrm{d}t$) shows a significant underestimation (${\text{fCnetProduct}}^{\ast }<\mathrm{d}\text{cProduct}/\mathrm{d}t$) for CESM2-WACCM and MIROC-ES2L. In CESM2-WACCM, this underestimation arises because the influxes (fDeforestToProduct and fHarvestToProduct) are much smaller than the outflux (fProductDecomp) – for example, 0.03 and 0.84 Gt C yr⁻¹ versus 3.6 Gt C yr⁻¹ during ssp585, on average. For MIROC-ES2L, the mismatch results from a non-zero fProductDecomp, while no fDeforestToProduct or fHarvestToProduct fluxes are reported (Fig. A1). ACCESS-ESM1-5, CMCC-CM2-SR5, CMCC-ESM2, MPI-ESM1-2-LR, and UKESM1-0-LL show nearly identical fCnetProduct* and $\mathrm{d}\text{cProduct}/\mathrm{d}t$ in all experiments. Other ESMs do not report product carbon pool sizes but do report certain fluxes into or out of the product pool (Fig. A1), resulting in a non-zero net flux (Fig. A7).

Figure A8 Comparison of net nitrogen fluxes for organic and mineral pools (fNnetOrganic* and fNnetMineral*, reconstructed from flux data) and the change of nitrogen pool size over time for organic and mineral pools ($\mathrm{d}\text{nOrganic}*/\mathrm{d}t$ and $\mathrm{d}\text{nMineral}/\mathrm{d}t$, derived from pool size data) from CMIP6 ESMs across different scenarios.

Figure A9 The mismatch of changes in land nitrogen pool size derived from pool size data (land nitrogen: nLand*(t) − nLand*(t0); organic nitrogen: nOrganic*(t) − nOrganic*(t0); mineral nitrogen: nMineral(t) − nMineral(t0)) (solid lines) and those retrieved from the integration of net nitrogen flux reconstructed from flux data (land nitrogen: fNnetLand*, organic nitrogen: fNnetOrganic*, and mineral nitrogen: fNnetMineral*, dashed lines) in CMIP6 ESMs across different scenarios. Note there is no nitrogen cycle in CanESM.

Figure A10 The carbon mass imbalance from mismatch of net land carbon flux (fCnetLand*, with or without the consideration of fLuc, reconstructed from flux data) and the change in land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) for CMIP6 ESMs across different scenarios.

Figure A11 The carbon mass imbalance from mismatch of net land carbon flux (fCnetLand*, with or without the consideration of fFire, reconstructed from flux data) and the change in land carbon pool size over time ($\mathrm{d}{\text{cLand}}^{\ast }/\mathrm{d}t$, derived from pool size data) for CMIP6 ESMs across different scenarios.