Synthesizing global carbon–nitrogen coupling effects – the MAGICC coupled carbon–nitrogen cycle model v1.0

The integration of a nitrogen cycle represents a recent advancement in Earth system models (ESMs). However, diverse formulations introduce uncertainty in the nitrogen effect on the carbon cycle, leaving the global carbon–nitrogen coupling effect unclear. In this study, we present CNit v1.0, a newly developed carbon–nitrogen cycle model designed for integration with MAGICC (Model for the Assessment of Greenhouse Gas Induced Climate Change), a widely used reduced-complexity model. CNit v1.0 has been calibrated to two land surface models (CABLE and OCN) and (the land component of) a set of Coupled Model Intercomparison Project Phase 6 (CMIP6) ESMs. CNit v1.0 is able to capture the dynamics of the more complex models' carbon–nitrogen cycle at the global-mean, annual scale. The emulation results suggest a consistent nitrogen limitation on net primary production (NPP) in CMIP6 ESMs, persisting throughout the simulations (i.e., over the period 1850–2100) in most models. The emulation provides a way to disentangle diverse nitrogen effects on carbon pool turnovers in CMIP6 ESMs, with our results suggesting that nitrogen deficiency generally inhibits litter production and decomposition while enhancing soil respiration (from a multi-model mean perspective). However, this disentanglement is limited due to a lack of simulations from CMIP6 ESMs which would allow us to fully separate the nitrogen and carbon responses. The results imply a potential reduction in land carbon sequestration in the future due to nitrogen deficiency. Future studies will use CNit to further investigate the carbon–nitrogen coupling effect, including uncertainty, in future climate projections.

1 Introduction

Atmosphere–ocean general circulation models (AOGCMs) and Earth system models (ESMs) are currently the most powerful tools we have for integrating our understanding of climate physics and providing comprehensive projections of the global climate and its variability (Meehl, 1990). However, these complex models require large computational power for their simulations, while the difference in assumptions, parameterizations, and structures across models often hinders a systematic quantification of uncertainties (Ohgaito et al., 2013). To combine the latest insights from various AOGCMs and ESMs, simple climate models (SCMs) – also called reduced-complexity models (RCMs) – are developed and routinely updated to represent and integrate the full uncertainty spectrum across the cause–effect chain of climate change (Nicholls et al., 2020, 2021). The highly parameterized formulations used in RCMs can, in some cases, parameterize the structural uncertainties from complex models. As a result, this flexibility of RCM structures allows for factor separation analysis to disentangle the key processes affecting the climate. With these features, RCMs are widely used for ensemble projections of scenarios and regularly feed into climate policy.

The Model for the Assessment of Greenhouse Gas Induced Climate Change (MAGICC), originally introduced by Wigley and Raper (1987, 1992, 2001) and further developed since (Meinshausen et al., 2011a, 2020; Nauels et al., 2017), is a key RCM that has been used for scenario classification in multiple IPCC reports (e.g., Climate Change 2014: Synthesis Report, IPCC, 2014; Global Warming of 1.5 °C, IPCC, 2018; Climate Change 2023: Synthesis Report, IPCC, 2023). MAGICC's main design principle is this: be as simple as possible while as mechanistic as necessary in the sense of being based on physical principles and/or long-term ESM calibrations (Meinshausen et al., 2011a).

The nitrogen cycle is a critical part in the Earth system's biogeochemistry and has a significant impact on climate alongside other element cycles like carbon or phosphorus (Fowler et al., 2013; Elser et al., 2007). As an essential nutrient for numerous fundamental biological processes, nitrogen is one of the major factors controlling the terrestrial carbon cycle and thus influences the carbon–concentration and carbon–climate feedbacks (Zaehle et al., 2010; Zaehle and Dalmonech, 2011; Fowler et al., 2013; Zaehle, 2013), the two main carbon cycle feedbacks (Arora et al., 2020). The integration of the nitrogen cycle and its effects within carbon cycle models is a recent advancement in ESMs. Only three Coupled Model Intercomparison Project Phase 5 (CMIP5) ESMs (CCSM, CESM, NorESM), all of which had the same land component (CLM4), included the nitrogen cycle (Flato et al., 2014). However, at least 17 out of 39 CMIP6 ESMs included a nitrogen cycle (see “Annex II: Models” in Climate Change 2021: The Physical Science Basis; IPCC, 2021). Various assumptions and formulations have been incorporated into the nitrogen cycle and the carbon–nitrogen coupling (Meyerholt and Zaehle, 2015; Meyerholt et al., 2020), resulting in divergent responses of the carbon cycle (Zaehle et al., 2015; Davies-Barnard et al., 2020; Arora et al., 2020; Kou-Giesbrecht and Arora, 2022).

Generally, the inclusion of the nitrogen cycle reduces land carbon sequestration under increasing atmospheric CO₂ and warming conditions by limiting photosynthesis (thus limiting net primary production, NPP) and amplifying both plant respiration and soil organic matter decomposition (Thornton et al., 2007; Sokolov et al., 2008). On average, the carbon–nitrogen coupled ESMs have smaller carbon–concentration feedbacks and smaller carbon–climate feedbacks (weaker absolute strength of the feedback parameters) compared to their carbon-only counterparts (Arora et al., 2020). Plant nitrogen uptake, the carbon : nitrogen ratio, the nitrogen regulation of photosynthesis, and biological nitrogen fixation contribute to the NPP difference (Du et al., 2018). Carbon–nitrogen interaction simulations from the Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg (JSBACH) have suggested a moderate reduction in the carbon–concentration feedback while showing a negligible effect on the carbon–climate feedback (Goll et al., 2017). However, enhanced soil organic matter decomposition under warming increases mineral nitrogen availability, thereby leading to increased land carbon sequestration on vegetation. The relative strength of these compensating effects remains unclear. Therefore, there is a need for better understanding and comparison of ESMs, with the integration and parameterization of a nitrogen cycle in RCMs providing one way to develop this understanding and comparison.

The significance of the nitrogen cycle also highlights the need to capture its effects within a key tool in climate science, namely RCMs. To the best of our knowledge, there is currently no RCM featuring the nitrogen effect, let alone a fully coupled carbon–nitrogen cycle. As a result, this study introduces a coupled carbon–nitrogen model (referred to as CNit), which has evolved from the previous MAGICC carbon cycle model. Section 2 presents a detailed description of the CNit model. Section 3 provides the results of offline calibration of CNit, firstly to two land surface models and then to a series of CMIP6 ESMs across multiple scenarios. Section 4 offers related discussions and analysis of the carbon–nitrogen coupling effect, primarily focusing on the CMIP6 ESMs. Section 5 discusses the limitations and implications of CNit's emulation. The results demonstrate that CNit captures the global aggregate effects of coupling the nitrogen cycle with the carbon cycle, in line with the latest generation of specialized domain models and ESMs. Future work will use MAGICC, updated to include CNit, to explore one of the key uncertainties in future climate projections: the uncertain evolution of future CO₂ concentrations given the intertwined carbon cycle feedback, CO₂ fertilization, and nitrogen cycle effects.

2 Model description

2.1 Overview of MAGICC and CNit

MAGICC is one of the most widely used RCMs. MAGICC features variable climate sensitivities and a carbon cycle that has successfully emulated a series of CMIP3 AOGCMs and Coupled Climate–Carbon Cycle Model Intercomparison Project (C4MIP) carbon cycle models (Meinshausen et al., 2011a). The most recent updates of MAGICC include the introduction of variable climate sensitivities and the updated carbon cycle (Meinshausen et al., 2011a); the incorporation of a sea level model (Nauels et al., 2017); and various improvements over time with regard to, for example, radiative forcing schemes (Meinshausen et al., 2020). In its latest calibration, MAGICC was shown to reproduce the IPCC AR6 WG1 assessment well over a range of metrics (see Chap. 7 in Climate Change 2021: The Physical Science Basis; IPCC, 2021).

Figure 1 The coupled carbon–nitrogen cycle model (CNit v1.0) in MAGICC (NPP: net primary production; LPR: litter production respiration; PU: nitrogen plant uptake; BNF: biological nitrogen fixation; LP: litter production; LD: litter decomposition; SR: soil respiration; LS: mineral nitrogen loss; LU: land use emission; AD: atmosphere nitrogen deposition; FT: nitrogen fertilizer application; 2P, 2L, 2S, and 2M: the partition of fluxes into plant, litter, soil, and mineral pools). The carbon and nitrogen mass conservation is described in Sect. 2.2. The NPP simulation is described in Sect. 2.3. The carbon–nitrogen coupling is described in Sect. 2.4. The LPR flux is described in Sect. 2.5. Each pool's turnover time and its response to climate and carbon–nitrogen coupling are described in Sect. 2.6. The land use emission parameterization is described in Sect. 2.7.

The continuously expanding understanding of climate physics, chemistry, and biology, coupled with the rapid development of complex models, necessitates corresponding advancement of RCMs. Here we focus on the development of the terrestrial carbon–nitrogen cycle model in MAGICC, referred to as CNit, which builds upon the previous MAGICC carbon cycle (Meinshausen et al., 2011a). As a brief background, the initial design of CNit considered carbon and nitrogen processes at a level of detail similar to that of complex models (albeit at a global scale rather than grid-box scale; e.g., the box model design starts from the major state variables and fluxes that are required by C4MIP) (Jones et al., 2016). However, during model parameterization and refinement, some processes were simplified or integrated with others to improve efficiency. For instance, biological nitrogen fixation is directly allocated to organic nitrogen pools, bypassing the intermediate step of mineral nitrogen enrichment and subsequent plant uptake (Fig. 1). Additionally, certain representations, such as land use emissions, were updated to achieve a balance between model simplicity and mechanistic insights (Sect. 2.7). These refinements align with MAGICC's need for computational efficiency and overall design philosophy of being as simple as possible but as mechanistic as necessary.

CNit is a globally integrated, annually averaged box model (Fig. 1) designed to simulate terrestrial carbon and nitrogen dynamics. It includes carbon and nitrogen pools for plant (P), litter (L), and soil (S), along with an inorganic mineral (M) nitrogen pool. The atmosphere (A) exchanges carbon with the land carbon pools via net primary production (NPP), heterotrophic respiration (RH), and land use or other anthropogenic fluxes (LU_c). Similarly, the atmosphere exchanges nitrogen with the land nitrogen pools via nitrogen atmospheric deposition (AD), biological nitrogen fixation (BNF), gaseous nitrogen loss (LS2A), and land use or other anthropogenic fluxes (LU_n). CNit takes the land use emissions of carbon and nitrogen, nitrogen fertilizer application (FT), AD, and BNF, as the inputs. Then it models key fluxes and solves a system of mass conservation equations to determine the fluxes and states for carbon and nitrogen. The resulting net land-to-atmosphere carbon and nitrogen fluxes are then used to estimate atmospheric concentrations, which subsequently inform radiative forcing and climate responses. These climate responses, in turn, interact with the carbon–nitrogen cycle, creating a feedback loop (see details in Meinshausen et al., 2011a).

The following sections outline the key components of CNit: Sect. 2.2 introduces the mass conservation framework and key fluxes; Sect. 2.3 details the NPP simulation; Sect. 2.4 explains carbon–nitrogen coupling, where we link nitrogen plant uptake (PU) and NPP; Sect. 2.5 describes the litter production respiration flux; Sect. 2.6 focuses on carbon and nitrogen turnover calculations; and Sect. 2.7 addresses the implementation of land use emissions.

2.2 Carbon and nitrogen mass conservation in CNit

The pools are interlinked by a system of first-order differential equations (Eqs. 1–9). In the equations, the subscripts “c” and “n” for the turnover fluxes (litter production (LP), litter decomposition (LD), and soil respiration (SR)) and land use fluxes (LU) denote the carbon and nitrogen fluxes. The subscripts X2Y for the partitioning factor “f” refer to the fraction of flux X that enters pool Y, whereas for land use flux, the subscripts LU2Y specifically represent the partitioning of the flux leaving pool Y due to land use changes. The flux partition accounts for the intended time domain of applicability – specifically, the long time step that integrates carbon and nitrogen cycle dynamics over periods of 1 year or more. Within this framework, fluxes such as NPP, BNF, PU, LU, and turnover fluxes not only contribute to changes in their immediate target pools but also propagate to subsequent pools. For instance, the NPP flux is simultaneously partitioned into plant, litter, and soil carbon pools (Eqs. 1–3). The partitioning factors always sum to unity, ensuring that no mass is lost as a result of this partitioning.

The carbon mass balance in plant (C_P), litter (C_L), and soil (C_S) pools is calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{P}}}{\mathrm{d}t}}=f_{\mathrm{NPP}\mathrm{2}\mathrm{P}}\mathrm{NPP}-\mathrm{LPR}-{\mathrm{LP}}_{\mathrm{c}}-f_{\mathrm{LU}\mathrm{2}{\mathrm{P}}_{\mathrm{c}}}{\mathrm{LU}}_{\mathrm{c}},\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{L}}}{\mathrm{d}t}}=f_{\mathrm{NPP}\mathrm{2}\mathrm{L}}\mathrm{NPP}+f_{\mathrm{LP}\mathrm{2}{\mathrm{L}}_{\mathrm{c}}}{\mathrm{LP}}_{\mathrm{c}}-{\mathrm{LD}}_{\mathrm{c}}-f_{\mathrm{LU}\mathrm{2}{\mathrm{L}}_{\mathrm{c}}}{\mathrm{LU}}_{\mathrm{c}},\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{S}}}{\mathrm{d}t}}& =f_{\mathrm{NPP}\mathrm{2}\mathrm{S}}\mathrm{NPP}+f_{\mathrm{LP}\mathrm{2}{\mathrm{S}}_{\mathrm{c}}}{\mathrm{LP}}_{\mathrm{c}}+f_{\mathrm{LD}\mathrm{2}{\mathrm{S}}_{\mathrm{c}}}{\mathrm{LD}}_{\mathrm{c}}-{\mathrm{SR}}_{\mathrm{c}}\\ & -f_{\mathrm{LU}\mathrm{2}{\mathrm{S}}_{\mathrm{c}}}{\mathrm{LU}}_{\mathrm{c}}.\end{array}\end{array}$$

Note that the litter production respiration (LPR) is the litter respiration produced from plant litter production that is released back to the atmosphere within a single time step (typically 1 year for MAGICC and CNit). Further details on the LPR flux are provided in Sect. 2.5.

For the total land carbon (sum of plant, litter, and soil carbon, i.e., combining Eqs. 1–3), we obtain

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{LAND}}}{\mathrm{d}t}}=\mathrm{NPP}-\mathrm{RH}-{\mathrm{LU}}_{\mathrm{c}}.\end{array}$$

Note that RH includes LPR, SR, and the litter decomposition that directly goes into atmosphere (LD2A, i.e., litter respiration).

The nitrogen mass balance in plant (N_P), litter (N_L), soil (N_S), and mineral (N_M) pools is calculated as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{P}}}{\mathrm{d}t}}=f_{\mathrm{BNF}\mathrm{2}\mathrm{P}}\mathrm{BNF}+f_{\mathrm{PU}\mathrm{2}\mathrm{P}}\mathrm{PU}-{\mathrm{LP}}_{\mathrm{n}}-f_{\mathrm{LU}\mathrm{2}{\mathrm{P}}_{\mathrm{n}}}{\mathrm{LU}}_{\mathrm{n}},\text{(6)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{L}}}{\mathrm{d}t}}& =f_{\mathrm{BNF}\mathrm{2}\mathrm{L}}\mathrm{BNF}+f_{\mathrm{PU}\mathrm{2}\mathrm{L}}\mathrm{PU}+f_{\mathrm{LP}\mathrm{2}{\mathrm{L}}_{\mathrm{n}}}{\mathrm{LP}}_{\mathrm{n}}-{\mathrm{LD}}_{\mathrm{n}}\\ & -f_{\mathrm{LU}\mathrm{2}{\mathrm{L}}_{\mathrm{n}}}{\mathrm{LU}}_{\mathrm{n}},\end{array}\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{S}}}{\mathrm{d}t}}& =f_{\mathrm{BNF}\mathrm{2}\mathrm{S}}\mathrm{BNF}+f_{\mathrm{PU}\mathrm{2}\mathrm{S}}\mathrm{PU}+f_{\mathrm{LP}\mathrm{2}{\mathrm{S}}_{\mathrm{n}}}{\mathrm{LP}}_{\mathrm{n}}\\ & +f_{\mathrm{LD}\mathrm{2}{\mathrm{S}}_{\mathrm{n}}}{\mathrm{LD}}_{\mathrm{n}}-{\mathrm{SR}}_{\mathrm{n}}-f_{\mathrm{LU}\mathrm{2}{\mathrm{L}}_{\mathrm{n}}}{\mathrm{LU}}_{\mathrm{n}},\end{array}\text{(8)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{M}}}{\mathrm{d}t}}=\mathrm{AD}+\mathrm{FT}+f_{\mathrm{LD}\mathrm{2}{\mathrm{M}}_{\mathrm{n}}}{\mathrm{LD}}_{\mathrm{n}}+{\mathrm{SR}}_{\mathrm{n}}-\mathrm{PU}-\mathrm{LS}.\end{array}$$

Note that (a) the mineral nitrogen loss (LS) is the mineral nitrogen turnover, (b) the nitrogen plant uptake (PU) is the mineral nitrogen taken up by the organic nitrogen pools, (c) the mineral nitrogen pool receives additional nitrogen from nitrogen fertilizer application (FT), and (d) the sum of the fraction of litter decomposition nitrogen entering the mineral pool (LD2M_n, i.e., litter mineralization) and the nitrogen released during soil respiration (SR_n, i.e., soil organic matter mineralization) constitutes the ecosystem's net mineralization (netMIN).

For the total land nitrogen (sum of plant, litter, soil, and mineral nitrogen, i.e., combining Eqs. 5–8), we have

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{LAND}}}{\mathrm{d}t}}=\mathrm{BNF}+\mathrm{AD}+\mathrm{FT}-\mathrm{LS}-{\mathrm{LU}}_{\mathrm{n}}.\end{array}$$

2.3 NPP simulation: effect of CO₂ and temperature forcings

The NPP flux is modeled by scaling an initial NPP (NPP₀) with the effect from changes in atmospheric CO₂ (${\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}$), temperature change (ϵ_dT(NPP)), carbon–nitrogen coupling (ϵ_CN(NPP)), and land use change (ϵ_LU):

$$\begin{array}{}\text{(10)}& \mathrm{NPP}={\mathrm{NPP}}_{\mathrm{0}}\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}\times {\mathit{ϵ}}_{\mathrm{CN}\left(\mathrm{NPP}\right)}\times {\mathit{ϵ}}_{\mathrm{LU}}.\end{array}$$

2.3.1 CO₂ fertilization

The CO₂ fertilization formulations can take multiple forms. The first, the logarithmic formulation, is adapted from Bacastow and Keeling (1973) as follows:

$$\begin{array}{}\text{(11)}& {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }=\mathrm{1}+s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }\times \ln \left({\mathrm{CO}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}\mathrm{ref}}\right),\end{array}$$

where $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }$ represents the sensitivity of NPP to the logarithm of the ratio of current atmospheric CO₂ concentration (CO₂) to a reference CO₂ level (CO_2ref, e.g., the pre-industrial CO₂ concentration) (i.e., the relative change from CO₂ to CO_2ref).

The second, the rectangular hyperbolic formulation, is adapted from Hunt et al. (1991) and Gifford (1993) as follows:

$$\begin{array}{}\text{(12)}& {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}={\displaystyle \frac{\mathrm{1}/\left({\mathrm{CO}}_{\mathrm{2}\mathrm{ref}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)+s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}}{\mathrm{1}/\left({\mathrm{CO}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)+s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}}},\end{array}$$

where CO_2b is the CO₂ concentration when NPP = 0, which has a default value of 31 ppm (Gifford, 1993), and $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}$ determines the CO₂ sensitivity of NPP in the rectangular hyperbolic formulation.

When the CO₂ concentration increases from 340 to 680 ppm, the ratio of the feedback factor at 680 ppm to that at 340 ppm (r) is designed to be the same for both formulations to ensure better compatibility:

$$\begin{array}{}\text{(13)}& r={\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }\left(\mathrm{680}\right)/{\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }\left(\mathrm{340}\right)={\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}\left(\mathrm{680}\right)/{\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}\left(\mathrm{340}\right).\end{array}$$

The sensitivities of NPP in the two formulations are therefore related by

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }r={\displaystyle \frac{\mathrm{1}+s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }\times \ln \left(\mathrm{680}/{\mathrm{CO}}_{\mathrm{2}\mathrm{ref}}\right)}{\mathrm{1}+s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }\times \ln \left(\mathrm{340}/{\mathrm{CO}}_{\mathrm{2}\mathrm{ref}}\right)}},\text{(15)}& {\displaystyle }& {\displaystyle }{s}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}={\displaystyle \frac{\left(\mathrm{680}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)-r\left(\mathrm{340}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)}{(r-\mathrm{1})\left(\mathrm{680}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)\left(\mathrm{340}-{\mathrm{CO}}_{\mathrm{2}\mathrm{b}}\right)}}.\end{array}$$

The previous MAGICC carbon cycle model uses a linear combination of the logarithmic and rectangular hyperbolic formulations to calculate the final CO₂ fertilization effect (Meinshausen et al., 2011a). However, because the logarithmic formulation is an unbounded function, the linear combination becomes unbounded as well unless the logarithmic formulation is removed, resulting in an overreliance on the rectangular hyperbolic formulation. The rectangular hyperbolic formulation itself can increase steeply if the CO₂ sensitivity of NPP in the rectangular formulation ($s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}$) is small. Considering $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}$ is dependent on r (or $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }$, Eqs. 14 and 15), the small $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }$ value is easily attainable, thereby leading to a high CO₂ fertilization factor. To fix this problem, a sigmoidal CO₂ fertilization formulation is introduced and included in the updated carbon–nitrogen model presented here.

$$\begin{array}{}\text{(16)}& {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sig}}={\displaystyle \frac{{\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}\max }^{\mathrm{sig}}}{\mathrm{1}+\left({\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}\max }^{\mathrm{sig}}-\mathrm{1}\right)\times e^{-s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sig}}\left({\mathrm{CO}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{ref}}\right)}}},\end{array}$$

where ${\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}\max }^{\mathrm{sig}}$ denotes the maximum of the sigmoidal CO₂ fertilization (always ≥1, which occurs when CO₂ reaches infinity) and $s_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sig}}$ is the CO₂ sensitivity of NPP in the sigmoidal formulation.

We allow for a linear combination of CO₂ fertilization formulations. A method factor ($m_{{\mathrm{CO}}_{\mathrm{2}}}$) ranging from 0 to 2 is used to combine the formulations and calculate the effective CO₂ feedback (illustrated in Fig. A1).

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}=\\ & \left\{\begin{array}{ll}\left(\mathrm{1}-m_{{\mathrm{CO}}_{\mathrm{2}}}\right)\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }+m_{{\mathrm{CO}}_{\mathrm{2}}}\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}& \mathrm{0}\le m_{{\mathrm{CO}}_{\mathrm{2}}}\le \mathrm{1}\\ \left(\mathrm{2}-m_{{\mathrm{CO}}_{\mathrm{2}}}\right)\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{rect}}+\left(m_{{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{1}\right)\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sig}}& \mathrm{1}<m_{{\mathrm{CO}}_{\mathrm{2}}}\le \mathrm{2}\end{array}\right.\end{array}\end{array}$$

2.3.2 Feedback from temperature change

Global-mean temperature change (dT) is taken as a proxy for climate-related impacts on the carbon cycle fluxes, i.e., for representing the carbon–climate feedback. The response of NPP to dT is assumed to follow an exponential or sigmoidal scaling of NPP₀, determined by a given temperature sensitivity ($s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\exp }$ or $s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\mathrm{sig}}$). The latter is introduced to better capture the trend of NPP in low-emission scenarios.

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\exp }=e^{s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\exp }\times \mathrm{d}T}\text{(19)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\mathrm{sig}}={\displaystyle \frac{\mathrm{2}}{\mathrm{1}+e^{-s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\mathrm{sig}}\times \mathrm{d}T}}}\end{array}$$

Similarly, a method factor (m_dT, between 0 and 1) is used to control the effective temperature change feedback on NPP (illustrated in Fig. A1).

$$\begin{array}{}\text{(20)}& {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}=\left(\mathrm{1}-m_{\mathrm{d}T}\right)\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\exp }+m_{\mathrm{d}T}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\mathrm{sig}}\end{array}$$

2.4 The carbon–nitrogen coupling in CNit

The current carbon–nitrogen coupling is based on the mineral nitrogen requirement and availability, from which the nitrogen deficiency (or surplus) can be calculated and, subsequently, the influence on NPP can be determined. A direct link between NPP (or photosynthesis) and plant nitrogen status is a common treatment in complex carbon–nitrogen coupled models (Zaehle and Dalmonech, 2011; Zaehle et al., 2014).

The overall formulation design of the carbon–nitrogen coupling effect on NPP is as follows: first, we establish a relationship between NPP (carbon fixation) and nitrogen plant uptake (PU, nitrogen fixation). We calculate the potential NPP (NPP_pot, i.e., the carbon-only NPP) by setting ϵ_CN(NPP)=1 in Eq. (10). The corresponding PU requirement (PU_req) is determined by this NPP_pot (as well as a temperature effect; see details in Eqs. 21 and 22). Then the nitrogen deficiency or surplus is calculated based on PU_req and nitrogen atmospheric deposition (AD), with which the ϵ_CN(NPP) is updated and the actual NPP (NPP_act) is determined. The following describes the process in detail.

2.4.1 The nitrogen plant uptake (PU) and NPP

First, we assume that the PU is a function of the NPP (whether potential or actual) and scaled by the temperature response of PU (ϵ_dT(PU)):

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }\mathrm{PU}={\mathrm{PU}}_{\max }\times e^{-\frac{{\mathrm{NPP}}_{\mathrm{ref}}}{\mathrm{NPP}}}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{PU}\right)},\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{PU}\right)}=e^{s_{\mathrm{d}T\left(\mathrm{PU}\right)}\times \mathrm{d}T},\end{array}$$

where NPP_ref is a reference NPP used to normalize the real-time NPP (further explanation in the next paragraph), PU_max sets the upper bound of the PU without the temperature feedback, and s_dT(PU) is the temperature sensitivity of PU.

When NPP=NPP_ref, with dT=0, the PU is fixed to PU${}_{\max }/e$. The PU_max and NPP_ref parameter pair defines a unique plant carbon–nitrogen assimilation system. Specifically, PU${}_{\max }/(e\times {\mathrm{NPP}}_{\mathrm{ref}})$ reflects a default setting of the plant nitrogen : carbon ratio. It is also worth reiterating that NPP_ref is not necessarily the same as the initial NPP (NPP₀ in Eq. 10).

The formulation presented in Eq. (21) suggests an increasing PU with increasing NPP. However, as NPP increases, the PU needed per additional unit of NPP gradually decreases. This provides a way for the model to represent a declining carbon : nitrogen ratio in plant biomass with increasing NPP due to, for example, CO₂ fertilization.

2.4.2 The potential NPP (NPP_pot) and plant nitrogen uptake requirement (PU_req)

When the nitrogen effect is not considered, NPP can reach NPP_pot, which can be calculated using Eq. (10) (fixing ϵ_CN(NPP)=1). The PU_req can be calculated using the universal PU formula (Eq. 21) by replacing NPP with NPP_pot:

$$\begin{array}{}\text{(23)}& {\mathrm{PU}}_{\mathrm{req}}={\mathrm{PU}}_{\max }\times e^{-\frac{{\mathrm{NPP}}_{\mathrm{ref}}}{{\mathrm{NPP}}_{\mathrm{pot}}}}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{PU}\right)}.\end{array}$$

The integration of PU_req over a certain time period (e.g., the model time step) indicates the amount of mineral nitrogen needed to realize the potential NPP.

2.4.3 The nitrogen effect on NPP

The mineral nitrogen availability depends on the current mineral nitrogen pool size and fluxes (Eq. 8). The net mineralization (netMIN), which is the largest natural source of mineral nitrogen, comes from litter and soil nitrogen turnovers. The turnovers are dependent on the pool sizes based on a first-order decay formulation (see Sect. 2.6). Considering that PU is the predominant influx for the organic nitrogen pools (5 to >10 times greater than biological nitrogen fixation in the complex models we have examined), the accumulation of nitrogen in plant, litter, and soil pools is effectively controlled by PU. In other words, PU and netMIN, the two major fluxes channeled through the organic and inorganic nitrogen pools either by consuming mineral nitrogen to enrich organic nitrogen or vice versa, are closely intertwined and mutually influence each other. This is supported by the results from complex models (CABLE, OCN, and multiple CMIP6 models) that show a similar value and trend for PU and netMIN at the global-mean, annual-mean level (Fig. A2).

When approximating netMIN as being linearly correlated with PU, the unmet PU_req from netMIN alone can be calculated as a linear function of the PU_req itself (i.e., $\mathrm{netMIN}-{\mathrm{PU}}_{\mathrm{req}}=f\left({\mathrm{PU}}_{\mathrm{req}}\right)$). Considering that nitrogen atmospheric deposition (AD) also supports the PU_req, its effect of alleviating the nitrogen deficiency is added by another linear function. The carbon–nitrogen coupling effect on NPP (ϵ_CN(NPP)) is then determined by

$$\begin{array}{}\text{(24)}& {\mathit{ϵ}}_{\mathrm{CN}\left(\mathrm{NPP}\right)}={\mathit{ϵ}}_{\mathrm{CN}\left(\mathrm{NPP}\right)\mathrm{0}}\times e^{f_{\mathrm{1}}\times \mathrm{AD}-f_{\mathrm{2}}\times {\mathrm{PU}}_{\mathrm{req}}},\end{array}$$

where ϵ_CN(NPP)0 is a base nitrogen effect on NPP when there is neither deficiency nor surplus. f₁ and f₂ are fitted parameters whose values are always positive. The ($f_{\mathrm{1}}\times \mathrm{AD}-f_{\mathrm{2}}\times {\mathrm{PU}}_{\mathrm{req}}$) term determines the relative strength of the current AD and unmet PU_req (i.e., nitrogen deficiency/surplus).

This formulation is transformed from complex models with the key idea of comparing mineral nitrogen availability and plant nitrogen requirement. In complex carbon–nitrogen models, the nitrogen availability is typically based on the current mineral nitrogen pool size (with mass unit), and the nitrogen requirement is computed from the integrated fluxes in a given time step (with mass unit) (Thornton et al., 2007; Wiltshire et al., 2021; Zaehle et al., 2014). The competition from microbial immobilization is also considered in some complex models. However, in a model with a much longer time step (e.g., annual) like ours, such a system would be inherently unstable since the mineral nitrogen pool size would be orders of magnitude smaller than the annual nitrogen demand (i.e., the system would be unstable because the turnover of the mineral nitrogen pool would be substantially smaller than the time step).

Fixing f₁=1 indicates that on the timescale of interest (e.g., annually), all mineral nitrogen from AD is 100 % available for the plant. However, this is not necessarily correct considering that (1) the process-level nitrogen limitation/fertilization does not remain constant over the timescale of interest (e.g., annually); (2) the mineral nitrogen accumulation from previous time steps can be used for the current time step; and (3) the direct mineral nitrogen loss, whose magnitude is determined by the mineral pool turnover time, may counterbalance the effect of nitrogen deposition on the fertilization. Therefore, giving freedom to both the f₁ and the f₂ parameters implicitly allows our formulation to consider the above effects.

Note that when both f₁ and f₂ parameters are calibrated, the ϵ_CN(NPP) does not necessarily have to be ϵ_CN(NPP)0 at the start of the experiments, which gives flexibility to the model to determine the carbon–nitrogen coupling effect in the pre-industrial condition.

2.4.4 The actual NPP and PU

After the nitrogen effect is calculated, the actual NPP (NPP_act, i.e., NPP with the carbon–nitrogen coupling effect) is determined by

$$\begin{array}{}\text{(25)}& {\mathrm{NPP}}_{\mathrm{act}}={\mathrm{NPP}}_{\mathrm{pot}}\times {\mathit{ϵ}}_{\mathrm{CN}\left(\mathrm{NPP}\right)}.\end{array}$$

And based on Eq. (21), the corresponding actual PU (PU_act) becomes

$$\begin{array}{}\text{(26)}& {\mathrm{PU}}_{\mathrm{act}}={\mathrm{PU}}_{\max }\times e^{-\frac{{\mathrm{NPP}}_{\mathrm{ref}}}{{\mathrm{NPP}}_{\mathrm{act}}}}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{PU}\right)}.\end{array}$$

The modeled NPP_act and PU_act are then used for the differential equations (Eqs. 1–9) to solve the carbon and nitrogen states.

2.5 The litter production respiration flux

MAGICC separately simulates a litter production respiration (LPR) flux (the fast litter respiration produced from the plant litter production that does not carry over into the subsequent year, i.e., that returns to the atmosphere on sub-annual timescales) by scaling an initial plant litter production respiration flux (LPR₀) with the CO₂ fertilization effect (${\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}$), the carbon–nitrogen coupling effect on NPP (ϵ_CN(NPP)), the climate effect (ϵ_dT(LPR)), and the land use change effect (ϵ_LU).

$$\begin{array}{}\text{(27)}& {\displaystyle }& {\displaystyle }\mathrm{LPR}={\mathrm{LPR}}_{\mathrm{0}}\times {\mathit{ϵ}}_{{\mathrm{CO}}_{\mathrm{2}}}\times {\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{LPR}\right)}\times {\mathit{ϵ}}_{\mathrm{CN}\left(\mathrm{NPP}\right)}\times {\mathit{ϵ}}_{\mathrm{LU}}\text{(28)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{d}T\left(\mathrm{LPR}\right)}=e^{s_{\mathrm{d}T\left(\mathrm{LPR}\right)}\times \mathrm{d}T}\end{array}$$

The plant litter production respiration flux is assumed to be a fast over-turning (e.g., within 1 year) outflux from the plant carbon pool that circulates through the NPP to plant to litter to the atmosphere within a single time step. Considering the close relationship between plant litter production respiration flux and NPP, it is scaled by the same CO₂ and carbon–nitrogen coupling feedback as NPP, as well as the same land use effect. Given that LPR results from three cascading processes – NPP, litter production, and decomposition – it is modeled with an exponential temperature response based on its own temperature sensitivity (s_dT(LPR)).

2.6 The turnover of carbon and nitrogen pools

The litter production (LP), litter decomposition (LD), and soil respiration (SR) for plant, litter, and soil carbon/nitrogen pools are assumed to be proportional to the corresponding pool sizes, linked by the turnover time (τ) and scaled by the effect from temperature change (ϵ_dT) and carbon–nitrogen coupling (ϵ_CN).

$$\begin{array}{}\text{(29)}& \mathrm{Turnover}={\displaystyle \frac{\mathrm{pool}{\mathrm{size}}_i}{{\mathit{\tau }}_i}}\times {\mathit{ϵ}}_{\mathrm{d}T}\times {\mathit{ϵ}}_{\mathrm{CN}},\end{array}$$

where i refers to the plant, litter, and soil carbon/nitrogen pools.

For the temperature feedback, it is assumed that each process has its own temperature sensitivity. The feedback is then modeled by an exponential relationship.

$$\begin{array}{}\text{(30)}& {\mathit{ϵ}}_{\mathrm{d}T\left(i\right)}=e^{s_{\mathrm{d}T\left(i\right)}\times \mathrm{d}T},\end{array}$$

where i refers to the carbon or nitrogen turnover process LP, LD, or SR.

The carbon–nitrogen coupling feedback takes current nitrogen plant uptake (PU) and nitrogen atmospheric deposition (AD) as proxies to represent the plant nitrogen status and the nitrogen forcing, respectively. And similarly, each turnover process has its own response to the current PU (s_PU(i)) and AD (s_AD(i)). An exponential relationship is used to simulate their effects on the processes.

$$\begin{array}{}\text{(31)}& {\mathit{ϵ}}_{\mathrm{CN}\left(i\right)}=e^{s_{\mathrm{PU}\left(i\right)}\times \mathrm{PU}}\times e^{s_{\mathrm{AD}\left(i\right)}\times \mathrm{AD}},\end{array}$$

where i refers to the carbon or nitrogen turnover processes LP, LD, or SR.

The mineral nitrogen loss (LS) is the turnover flux for the mineral nitrogen pool, which follows a first-order decay formulation (Eq. 29). However, only a temperature effect is applied to LS since the current AD and PU, the two proxies for the carbon–nitrogen feedback, are already directly the influx and outflux for the mineral pool (Fig. 1).

2.7 The updated implementation of land use emissions and their impact on NPP

Most scenarios run in MAGICC, for example those in AR6 (see Chap. 7 in Climate Change 2022: Mitigation of Climate Change; Riahi et al., 2022), directly report net land use emissions (gross deforestation − regrowth) instead of the separation of those two parts. Considering a box model with fixed turnover times, when the input is constant (e.g., constant NPP), a one-off subtraction of carbon out of one of the pools (e.g., one-off gross deforestation) would not lead to a permanent reduction in the amount of carbon in this pool. Instead, it would yield an asymptotic rebound to the pre-intervention carbon pool size over time (as the carbon that was moved out of the pool is taken up by the land again), which implies full regrowth in the long term. In MAGICC's previous carbon cycle, the partial regrowth was accounted for using a factor that determines the fraction of expected regrowth and adjustments to the turnover times (Meinshausen et al., 2011a). This implementation blends the regrowth with changes in the pool's turnover fluxes. As the turnovers are also impacted by feedbacks (Eq. 29), the previous MAGICC carbon cycle setup required a parallel calculation of the non-feedback case to retrieve the regrowth flux and correctly handle the land use input.

2.7.1 The updated land use emission implementation

The updated land use emission implementation described here makes the regrowth and gross deforestation more straightforward. The LU input in Eqs. (1)–(4) refers to the net land use emissions, which can be written as

$$\begin{array}{}\text{(32)}& {\mathrm{LU}}_{\mathrm{c}}={\mathrm{LU}}_{\mathrm{grsd}}-{\mathrm{LU}}_{\mathrm{rgr}},\end{array}$$

where LU_grsd and LU_rgr refer to the gross deforestation (instantaneous biomass extraction from land organic carbon pool) and regrowth (legacy biomass addition to land organic carbon pool), respectively.

The regrowth formulation assumes a constant regrowth flux during the growing years following each instance of gross deforestation. Therefore, the regrowth flux is calculated as the total regrowth (part of the gross deforestation) divided by the regrowth time, formulated as follows:

$$\begin{array}{}\text{(33)}& {\mathrm{LU}}_{\mathrm{rgr}}={\displaystyle \frac{\mathit{\phi }\underset{\max \left(\mathrm{0},t-{\mathit{\tau }}_{\mathrm{rgr}}\right)}{\overset{t}{\int }}{\mathrm{LU}}_{\mathrm{grsd}}}{{\mathit{\tau }}_{\mathrm{rgr}}}},\end{array}$$

where φ refers to the fraction of gross deforestation that can regrow and τ_rgr denotes the regrowth time required to reach the partial regrowth.

The formulation assumes a constant regrowth flux (LU_rgr) for τ_rgr years after every single gross deforestation occurrence – the regrowth flux at time t is thus the potential total regrowth (φ times the integration of LU_grsd) divided by the required regrowth time (τ_rgr).

Note that, in this formulation for regrowth and gross deforestation, both fluxes exist only within a certain time period (τ_rgr); hence they cannot change the long-term equilibrium of the system because the fixed turnover times always return the system to its pre-deforestation state (Eqs. 1–4). We discuss our solution for this in the next section.

2.7.2 The effect of land use emission on NPP

Instead of adjusting the turnover times (as in previous versions of the MAGICC carbon cycle), CNit applies the effect of land use change to NPP to change the long-term equilibrium after the gross deforestation (ϵ_LU in Eq. 10), which is formulated as

$$\begin{array}{}\text{(34)}& {\mathit{ϵ}}_{\mathrm{LU}}={\displaystyle \frac{C_{{\mathrm{LAND}}_{\mathrm{0}}}-(\mathrm{1}-\mathit{\phi })\underset{\mathrm{0}}{\overset{t}{\int }}{\mathrm{LU}}_{\mathrm{grsd}}}{C_{{\mathrm{LAND}}_{\mathrm{0}}}}},\end{array}$$

where $C_{{\mathrm{LAND}}_{\mathrm{0}}}$ is the initial land carbon pool size.

The land use effect on NPP assumes that the NPP is reduced immediately after gross deforestation activities. Considering the simultaneous regrowth, the NPP reduction is proportional to the accumulated permanent gross deforestation ($(\mathrm{1}-\mathit{\phi })\underset{\mathrm{0}}{\overset{t}{\int }}{\mathrm{LU}}_{\mathrm{grsd}}$). With the land use effect on NPP, when there is a one-off gross deforestation occurrence, the NPP gradually decreases to a new value, with that value depending on the regrowth fraction φ (e.g., φ=0 refers to zero regrowth and the gross deforestation causes permanent carbon loss from the land carbon pools). Because of the NPP change, the long-term equilibrium of the system changes too.

The new formulation enables the separate calculation of regrowth and gross deforestation without the non-feedback run. Arguably, it also makes more physical sense because deforestation reduces the amount of forest (area) available to grow and hence the amount of NPP. The regrowth amount is a linear function of regrowth time, representing a simplified approach compared to sigmoidal growth models. This is a trade-off that could be re-considered in future work. The applied effect on NPP also means an additional constraint on MAGICC's simulation of regrowth and gross deforestation. However, it should be noted that the definition discrepancy of land use change emissions itself still exists among different models/approaches and leads to substantial differences in land use emission estimates (Gasser and Ciais, 2013; Stocker and Joos, 2015; Grassi et al., 2023). Thus, the LU input must align with the outputs of complex models, the scenarios being applied, and specific definitions of land use change emissions (e.g., ELUC in the Global Carbon Budget project) to ensure consistency in model applications.

2.7.3 The limitation of land use emission in CNit

The formulation presented here has certain limitations. For example, it does not account for variations in regrowth rates among different ecosystems or across ecosystem successional stages – an inherent constraint of the global box model approach. Additionally, it aggregates deforestation and harvest fluxes into a single LU input, although the regrowth fraction parameter may provide some indication of harvest activities that do not result in regrowth.

It should also be noted that MAGICC, whether using the previous carbon cycle model or the CNit version, does not simulate carbon or nitrogen storage in the product pool due to its relatively small size (Eqs. 1–9 and Fig. 1). Thus, the deforestation and harvest fluxes going into the product pool, expected to be small, are accounted for in the land use emission input (LU_c or LU_n) and associated regrowth parameter, while the land carbon (or nitrogen) pool does not include the correspondingly stored carbon (or nitrogen) within the product.

3 Model calibration

3.1 Overview of the calibration process and results

This section presents the offline calibration results for CNit, i.e., calibrations using prescribed land surface temperature and atmospheric CO₂ concentration from the original model outputs. We first describe the data acquisition and post-processing of land surface model outputs and CMIP6 ESM outputs (Sect. 3.2). Next, we define the calibration targets (major fluxes and pool sizes) and weight them to create a cost function. Finally, we apply optimization algorithms to identify the “best-estimate” parameter set (Sect. 3.3). For a single model, all experiments are calibrated simultaneously, resulting in one best-estimate parameter set that captures the model's behavior across experiments. Using these best-estimate parameter values, we evaluate CNit emulation against model outputs and calculate the root mean square error (RMSE) and normalized RMSE to assess emulation performance. A discussion of the calibration results for CABLE, OCN, and CMIP6 ESMs is provided in Sect. 3.4 and 3.5.

3.2 Data acquisition and processing

3.2.1 Land surface models and outputs

The CABLE and OCN land surface model output datasets (global-mean, annual-mean values) were obtained directly from the modeling groups. CABLE is the land surface component for the Australian Community Climate and Earth System Simulator (ACCESS-ESM1) (Law et al., 2017) and ACCESS-CM2 (Bi et al., 2020). OCN is the updated land surface model built on ORCHIDEE (Zaehle and Friend, 2010), the land surface component of the IPSL Climate Model (IPSL-CM; Boucher et al., 2020). Both CABLE and OCN provided the results from the carbon-only and the carbon–nitrogen coupled setups for the Representative Concentration Pathway 8.5 (RCP8.5) scenario (Fleischer et al., 2019), with OCN also providing the results for the RCP2.6 scenario (Meyerholt et al., 2020). The data from both the carbon-only and carbon–nitrogen coupled setups of the land surface models offer valuable information to constrain nitrogen interactions separately from the climate and CO₂ effects. The carbon–nitrogen coupled runs in CABLE consisted of two experiments with constant and dynamic atmospheric deposition inputs, respectively (Fleischer et al., 2019). These experiments are useful for diagnosing the standalone effect of atmospheric deposition. The climate data and atmosphere CO₂ concentration for the land surface model experiments were derived from either their corresponding ESM outputs (Meyerholt et al., 2020) or the RCP greenhouse gas concentrations (Meinshausen et al., 2011b). One OCN model structure with a flexible carbon : nitrogen ratio, a linear biological nitrogen fixation and actual evaporation relationship, and explicit mineral nitrogen loss representation – namely the FLX/FOR/NL1 structure – from the 30 ensemble model structures in the original paper (Meyerholt et al., 2020) was selected for the CNit calibration in this paper. The selected OCN model structure serves as a proof of concept for the proposed CNit model. Future work will explore alternative structures and address structural uncertainties.

3.2.2 CMIP6 ESMs and outputs

We select target ESMs based on the following criteria: (1) the model includes a terrestrial nitrogen cycle; (2) the outputs encompass the majority of carbon–nitrogen cycle variables needed for C4MIP; and (3) the model has completed all selected experiments, including one idealized experiment (1pctCO2), historical simulations (historical), and four Shared Socioeconomic Pathway (SSP) scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5), all of which are CMIP6 tier 1 concentration-driven experiments. Based on the data availability and completeness (in the Earth System Grid Federation), six ESMs were chosen for this calibration (Table 1).

Table 1 The list of Earth system models (ESMs) selected for the calibration.

The monthly gridded CMIP6 ESM outputs were downloaded from the Earth System Grid Federation (ESGF, https://esgf-node.llnl.gov/projects/cmip6/, last access: 18 January 2024) and processed into global-mean, annual-mean values. Our grid-to-global aggregation used the model-specific grid area (areacella) and land fraction (sftlf) to avoid issues with the resolution variation across different model outputs. Our monthly-to-annual aggregation first concatenated the original monthly outputs along the time dimension and then calculated the annual values, weighted by the number of days in each month as defined by each ESM's output calendar. It should be noted that, even though the selected ESMs and experiments provided relatively complete outputs for the pools and fluxes, none of them reported all the required fluxes by C4MIP, especially those related to the land use and anthropogenic perturbations (Tang et al., 2025). As part of the processing, it was not trivial to reproduce the models' mass balance based on the reported outputs (Tang et al., 2025). For this paper, we have applied a workaround, which we describe in the next section.

3.3 Calibration setup

3.3.1 Temperature profile and CO₂ concentration

For CABLE and OCN, temperature and CO₂ concentration data were obtained directly from the original providers (Fig. A3). For CMIP6 ESMs, temperature data (tas) were sourced from model outputs. CNit uses temperature change (dT) as a proxy for climate-related impacts, calculated as the difference between the current temperature and the initial year's temperature (Figs. A3 and A4). CO₂ concentrations were taken from the CMIP6 forcing datasets (Meinshausen et al., 2017, 2020) – as all calibrated experiments are concentration-driven.

3.3.2 Inputs for CNit

The inputs for CNit (Fig. 1) include nitrogen atmospheric deposition, biological nitrogen fixation, and fertilizer application fluxes, all directly available from model outputs. For land use emissions of carbon and nitrogen, CABLE and OCN provide readily available data. However, CMIP6 ESMs do not provide such data (Tang et al., 2025). To address this, we calculated land use carbon emission (LU_c in Fig. 1 and Eqs. 1–4) as $\mathrm{NPP}-\mathrm{RH}-\mathrm{d}{C}_{\mathrm{LAND}}/\mathrm{d}t$, ensuring carbon mass conservation (as there is no straightforward land use emission from the reported outputs alone; see details in Tang et al., 2025).

Where available, land use and anthropogenic disturbance-related nitrogen fluxes (fNAnthDisturb and fNProduct in the CMIP6 Data Request) were combined to create the land use nitrogen input for CNit (LU_n in Fig. 1 and Eqs. 5–9).

3.3.3 Calibration target and optimization

Calibration targets for both land surface models and CMIP6 ESMs included NPP, heterotrophic respiration, nitrogen plant uptake (PU), and all carbon and nitrogen pool sizes. The cost function was calculated as the sum of normalized errors for each target flux or pool size time series (i.e., square $(\mathrm{emulation}-\mathrm{target})/({\mathrm{target}}_{\max }-{\mathrm{target}}_{\min })$). This normalization accounted for the differing magnitudes among target variables and gives lower weight to variables with greater variability. All available experiments were calibrated simultaneously without additional weighting, meaning the final cost was calculated as the sum of the costs from all experiments.

For CMIP6 ESMs, historical-period data and SSP scenario results were combined (referred to as hist_SSP) into a unified time axis spanning 1850–2100 or 1850–2300, depending on data availability.

The calibration process employed the differential evolution algorithm (Storn and Price, 1997) for global optimization and the Nelder–Mead algorithm (Gao and Han, 2012) for local minimization. First, differential evolution ran for 30 000 iterations with 10 random initializations. Next, the resulting 10 global-minimum parameter sets were used as initial guesses for the Nelder–Mead algorithm, yielding respective local minimums. Finally, the parameter set with the lowest cost was selected as the best-estimate parameter set. A full parameter list and the best-estimate values are provided in Tables A1 and A2.

3.4 Calibrating CNit to CABLE and OCN: results and comparison

CNit has successfully emulated the fluxes and pool sizes from the CABLE and OCN experiments (Figs. 2 and A5), with the RMSE (or normalized RMSE) ranging from 0.9–2.6 GtC yr⁻¹ (1.3 %–2.2 %) for NPP, 4.0–32.9 GtC (0.2 %–1.0 %) for the land carbon pool size, 0.016–0.020 GtN yr⁻¹ (1.1 %–1.9 %) for PU, 0.05–0.11 GtN (<0.1 %) for the land organic nitrogen pool size, and 0.0022–0.026 GtN (2.3 %–6.2 %) for the mineral nitrogen pool size. There is a strong and increasing nitrogen limitation on NPP in CABLE, inhibiting NPP by 9.4 % to 48.2 % (constant atmospheric deposition) or 46.2 % (dynamic atmospheric deposition) from 1901 to 2100. The total land carbon storage, i.e., the total amount of carbon taken up by the land carbon cycle, by 2100 under the RCP8.5 scenario in the constant and dynamic atmospheric deposition experiments is 353 and 398 GtC, respectively, which is significantly lower than the 1634 GtC land carbon storage in the runs that use a carbon cycle without nitrogen limitation. The large reduction in NPP (influx of the system, Fig. 2) due to the carbon–nitrogen coupling, along with the relatively small reduction in heterotrophic respiration (major outflux of the system, Fig. A5), results in the decrease in net land carbon flux (Fig. 2) and, thus, the decreased land carbon storage (Fig. 2).

Figure 2 Comparison of net primary production (NPP), land carbon pool size (C_LAND), net land carbon flux (fC_NetLAND), nitrogen plant uptake (PU), land organic nitrogen pool size (N_LANDo, sum of nitrogen in plant, litter, and soil pools), mineral nitrogen pool size (N_M), and net land nitrogen flux (fN_NetLAND) between CABLE or OCN outputs (blue lines) and CNit emulations (orange lines). The experiments labeled C, CN, and CNd denote the carbon-only, carbon–nitrogen coupled with constant nitrogen atmospheric deposition, and carbon–nitrogen coupled with dynamic nitrogen atmospheric deposition configurations, respectively, in the land surface models.

The OCN model exhibits significantly less nitrogen limitation on NPP (<5.0 % inhibition) for both the RCP8.5 and the RCP2.6 scenarios. This disparity may be attributed to the fact that the NPP in OCN in the carbon-only simulations (∼85 GtC yr⁻¹ in 2100 in RCP8.5) is notably lower than that in CABLE (∼140 GtC yr⁻¹ in 2100 in RCP8.5). The carbon cycle in OCN has also experienced a higher temperature change (∼8 °C from 1850 to 2100 in RCP8.5, Fig. A3) compared to CABLE (∼5.5 °C from 1900 to 2100 in RCP8.5, Fig. A3). The emulation captures the varied NPP increases by applying a stronger CO₂ fertilization effect and a more negative climate feedback in CABLE than in OCN. These respective changes are captured in CNit by the logarithmic carbon–concentration feedback formulation with a high sensitivity ($s_{{\mathrm{CO}}_{\mathrm{2}}}^{\log }=\mathrm{2.582}$) and negative $s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\exp }$ and $s_{\mathrm{d}T\left(\mathrm{NPP}\right)}^{\mathrm{sig}}$ parameters (Table A2). Nevertheless, the relatively minor NPP limitation in OCN, along with the resulting decrease in net land carbon flux (Fig. 2), leads to a reduction of 45 GtC (or 10 %, RCP8.5) and 26 GtC (or 8 %, RCP2.6) in land carbon storage by the end of 2100 compared to their respective simulations that use a carbon cycle without a nitrogen limitation.

In the RCP8.5 experiment, CABLE and OCN show similar relative changes in PU over time; however, their starting PU fluxes differ (∼0.6–1.0 GtN yr⁻¹ in CABLE vs. ∼1.0–1.4 GtN yr⁻¹ in OCN). The emulation is able to capture both dynamics by adjusting the maximum PU (PU_max) and temperature sensitivity of PU (s_dT(PU)) parameters. Specifically, the different starting PU values stem from the higher PU_max in OCN (2.4 GtN yr⁻¹) than in CABLE (1.9 GtN yr⁻¹). The similar trends occur because of the higher s_dT(PU) in CABLE compared to OCN (0.014 K⁻¹ vs. 0.008 K⁻¹, Table A2), which compensates for the effect of its lower PU_max (Eqs. 21 and 22). In the dynamic atmospheric deposition experiment, the PU in CABLE is slightly higher than in the constant atmospheric deposition experiment, primarily accumulating in the soil organic nitrogen pool (+3 GtN in 2100 compared to the constant atmospheric deposition scenario, Fig. A5). The accumulated organic nitrogen enriches the mineral nitrogen pool instead of facilitating the nitrogen loss, as the system remains nitrogen limited. In contrast, in the RCP8.5 experiment from OCN, the new organic nitrogen is mainly stored in the plant pool (Fig. A5). The trend and magnitude of the mineral nitrogen pool size in OCN are also largely different from those in CABLE. The emulation captures the diverse mineral nitrogen trends by adapting the temperature response of mineral nitrogen loss for the two models (weak and negative s_dT(LS) of −0.007 K⁻¹ in CABLE vs. strong and positive s_dT(LS) of 0.088 K⁻¹ for OCN, Table A2). The order-of-magnitude difference in mineral nitrogen pool sizes demonstrates the huge uncertainty in estimated mineral nitrogen quantities.

Comparing the behavior of the nitrogen cycle in OCN's RCP8.5 and RCP2.6 experiments, it is found that PU follows the trend of NPP, supporting our assumption that PU can be modeled based on NPP (Eq. 21). However, unlike in the RCP8.5 scenario, the new nitrogen introduced by PU is mainly accumulated in the soil nitrogen pool in the RCP2.6 scenario (Fig. A5). Based on the formulation, the larger soil nitrogen pool size implies more nitrogen mineralization (the first-order turnover, Eq. 29), enabling the emulation to capture the increasing trend of the mineral nitrogen pool size from 2050 to 2100 (Fig. A5). The organic nitrogen accumulation in the RCP8.5 and RCP2.6 experiments does not precisely follow the corresponding carbon storage trends, where the new carbon introduced by NPP is predominantly stored in the plant carbon pool in both scenarios. This difference indicates that the carbon cycle and nitrogen cycle, even though closely intertwined, may react in a divergent manner compared to climate change.

3.5 Calibrating CNit to CMIP6 ESMs: results and comparison

CNit emulation captures the dynamics of major fluxes and pool sizes for all the CMIP6 ESMs and experiments (Figs. 3 and A6), despite the diverse climate and CO₂ forcings (Fig. A4). More detailed comparisons of ESM outputs and CNit emulations are provided in Fig. A6. Overall, the model is able to emulate the wide range of ESM behavior. In the rest of this section, we focus on cases where this is not the case.

Figure 3 Comparison of net primary production (NPP), land carbon pool size (C_LAND), net land carbon flux (fC_NetLAND), nitrogen plant uptake (PU), land organic nitrogen pool size (N_LANDo, sum of nitrogen in plant, litter, and soil pools), mineral nitrogen pool size (N_M), and net land nitrogen flux (fN_NetLAND) between CMIP6 ESM outputs and CNit emulations. Results are normalized to a range of 0–1 using the following transformation: $x\mathrm{axis}=(\mathrm{target}-{\mathrm{target}}_{\min })/({\mathrm{target}}_{\max }-{\mathrm{target}}_{\min })$ and $y\mathrm{axis}=(\mathrm{emulation}-{\mathrm{target}}_{\min })/({\mathrm{target}}_{\max }-{\mathrm{target}}_{\min })$. The diagonal dashed line represents points where the emulation matches the target exactly, with positions below and above the line indicating underestimation and overestimation, respectively, by the emulator.

In the 1pctCO2 experiment, the emulated NPP is lower than that from the UKESM1-0-LL and MIROC-ES2L outputs. The underestimation is mainly from the inconsistent behavior of these two ESMs in the idealized 1pctCO2 experiment and hist_SSP experiments. Both ESMs have simulated higher NPP at the end of their 1pctCO2 runs (∼100 GtC yr⁻¹ in 1999 for both ESMs) than that at the end of their SSP runs (e.g., SSP1-2.6, <80 GtC yr⁻¹ in 2100 for both ESMs), which is contradictory to their PU results (lower in 1pctCO2 and higher in SSPs, Fig. A6). Such behavior is in direct contradiction with our assumption that higher NPP requires a higher PU (Sect. 2.4, Eq. 21). Since NPP and PU are both set as calibration targets, CNit has tried to minimize the gap between the emulated fluxes and targets, resulting in the simultaneous underestimation of NPP and overestimation of PU for UKESM1-0-LL and MIROC-ES2L (Figs. 3 and A6). However, such different behavior is only observed in the 1pctCO2 experiments for these two ESMs, indicating there could be either some model response nonlinearities between their 1pctCO2 and SSP runs that our model is not capturing or some regionally distinct effects that we are not seeing in the global, annual averages.

The underestimated NPP in the 1pctCO2 experiment has led to an underestimation of ∼190 GtC of the land carbon pool size for MIROC-ES2L (Fig. A6). But the underestimation of emulated land carbon in UKESM1-0-LL is much smaller (∼78 GtC, Fig. A6). This is because the CNit emulation has overestimated the soil carbon storage in the later phase of the 1pctCO2 experiment, which compensates for some carbon loss by the underestimated NPP. The emulated land carbon pool sizes (specifically the plant carbon pool sizes) are slightly but systematically smaller than the outputs from CMCC-CM2-SR5 and CMCC-ESM2 in their 1pctCO2 runs (Fig. A6). Such results primarily stem from CNit's underestimation of NPP during the middle of the 1pctCO2 experiment (Fig. A6). The CNit emulation of pool turnovers also contributes to the inconsistency between the emulation results and ESM outputs. The first-order decay of carbon pool turnover (Eq. 29) is not particularly effective in modeling the minor changes in pool sizes relative to the substantial initial pool size. For example, both CMCC models exhibit a soil carbon loss of approximately 45 GtC in their 1pctCO2 runs, with an initial soil pool size of around 2870 GtC. Consequently, the calibration has prioritized achieving a better fit for soil pool sizes at the expense of accurately representing plant and litter pool sizes. The different starting soil pool sizes for the ESMs' 1pctCO2 and historical simulations further complicate the soil carbon turnover emulation (Fig. A6). These different initial pool sizes in CMIP6 data do not make much sense. However, without explicit reasons to discredit the data, we have not adjusted our model to accommodate such oddities; rather, they compromise the model's fit.

The net land carbon flux is a highly variable flux (<10 GtC yr⁻¹ for all the ESMs, Fig. A6). CNit emulation of the ESMs has both positive and negative errors (Fig. 3). The CNit emulation has captured both their trends and magnitudes (Fig. A6).

The new nitrogen cycle in MAGICC has demonstrated the ability to capture the dynamics of nitrogen fluxes and pools, especially the PU and land organic nitrogen pool, in the ESM's SSP scenario runs (Fig. 3). In the 1pctCO2 experiment, however, the emulated PU is overestimated for UKESM1-0-LL and MIROC-ES2L, which is accompanied by the underestimation of their NPP, in an attempt to compensate for the conflicting high NPP and low PU in the ESMs' outputs (Fig. A6). The emulated organic nitrogen pool size is overestimated for NorESM2-LM (Fig. A6), which is primarily attributed to the overestimated PU and soil nitrogen pool sizes (Fig. A6).

The mineral nitrogen pool sizes exhibit the most diverse results between the CNit emulation and the ESM outputs (Fig. 3, details in Fig. A6). Firstly, the mineral nitrogen pool sizes in UKESM1-0-LL are relatively well emulated. Secondly, the emulated mineral nitrogen pool sizes show significant differences compared to NorESM2-LM outputs, in both magnitudes and trends. Considering that both the PU and the organic nitrogen pool size in NorESM2-LM are relatively well emulated (Fig. A6), such differences in the mineral nitrogen pool size could be attributed to either the input nitrogen fertilizer application flux or the turnover flux of the mineral pool. The CNit emulation of mineral nitrogen turnover is a first-order decay formulation with an exponential temperature effect scaling (Eq. 29). Thus, the initial pool size and turnover time determine the magnitude of the nitrogen loss flux. Among all the studied CMIP6 ESMs, NorESM2-LM has the largest initial mineral nitrogen pool (5.64 GtN vs. <1.50 GtN for all the other ESMs), which can lead to a large turnover flux. Conversely, the pool size change of mineral nitrogen is not significant in NorESM2-LM (0.78 GtN for 1pctCO2 and <0.35 GtN for all the hist_SSP experiments throughout the entire duration of the simulation). The large flux and small pool size change are naturally incompatible with each other in a first-order decay assumption, indicating uncaptured nonlinearities of the mineral pool simulation in NorESM2-LM. And lastly, the CNit emulation has effectively captured the trends in mineral nitrogen pool sizes simulated in CMCC-CM2-SR5, CMCC-ESM2, MPI-ESM1-2-LR, and MIROC-ES2L across all experiments (Fig. A6). Overall, the mineral pool is only a small part of the land nitrogen pool (the mineral-to-organic nitrogen ratio across the time series was 0.23 %–0.31 %, 0.24 %–0.31 %, 1.36 %–2.80 %, 2.05 %–2.53 %, 0.06 %–0.19 %, and 0.14 %–0.43 % for CMCC-CM2-SR5, CMCC-ESM2, MPI-ESM1-2-LR, NorESM2-LM, UKESM1-0-LL, and MIROC-ES2L, respectively). As a result, the net land nitrogen flux is mainly controlled by the organic nitrogen dynamics and shows consistency for the results from CNit emulation and ESM simulation (Figs. 3 and A6).

4 Discussions

4.1 The climate and carbon–nitrogen cycle in CMIP6 ESMs

The climate and carbon–nitrogen cycle in CMIP6 ESMs remain considerably different from each other (Figs. 4 and A7 and Sect. A1), which contributes to the imperfect emulation (especially the mineral nitrogen pool size, Fig. 3). The following discussion of the ESM outputs provides an overview of the current model-based understanding of the climate and carbon–nitrogen cycle from a global-mean, annual-mean perspective, from which we suggest future research needs. For the sake of comparison, the subsequent figures and discussions focus on the common experimental periods for scenarios and ESMs (e.g., 1850–2100 for the hist_SSP experiments). If not specified, the value and spread in the discussions are expressed as the mean ± 1 standard deviation across ESMs.

Figure 4 Carbon and nitrogen pool size dynamics (a: dC_LAND, delta land carbon pool size; b: dN_LANDo, delta land organic nitrogen pool size (sum of nitrogen in plant, litter, and soil pools); c: dN_M, delta mineral nitrogen pool size) from CMIP6 ESMs across different scenarios.

As a start, different ESMs have different temperature responses (Fig. A3), which has flow-on effects for the carbon cycle and nitrogen cycle (Fig. A6). Regardless of the absolute land carbon/nitrogen pool sizes (Fig. A8 and Sect. A2), the accumulation of land carbon/nitrogen during the same experimental period and their trends are considerably more similar (Fig. 4a and b). The accumulated land carbon storage is 472±201, 127±28, 164±50, 165±100, and 227±100 GtC for the 1pctCO2, hist_SSP126, hist_SSP245, hist_SSP370, and hist_SSP585 experiments, respectively. The corresponding land nitrogen accumulation is 3.04±4.65, 3.28±1.03, 3.92±1.33, 3.76±1.32, and 3.54±1.64 GtN. Generally, the land carbon accumulation and nitrogen accumulation are proportional to maintain the stoichiometric relationship between carbon and nitrogen (Fig. 4a and b, the mean values). However, the opposite trend of land carbon and nitrogen change is found in the two CMCC ESMs (compared to the other ESMs) in their 1pctCO2 runs, which explains the much larger spread of land organic nitrogen accumulation at the end of 1pctCO2 compared to the four hist_SSP experiments (Fig. 4b). Both the land carbon storage and the nitrogen storage show larger spread in the 1pctCO2 experiment (without land use changes), which further highlights the model structure uncertainty (and potentially also points to some inconsistency in how models are running these experiments). The higher spread of carbon and nitrogen storage in higher-warming scenarios indicates the feedback uncertainty (Melnikova et al., 2021).

The continuous and rapid depletion of mineral nitrogen is observed in NorESM2-LM under the 1pctCO2 scenario, coinciding with the highest accumulation of organic nitrogen (Fig. 4c). There are contrasting trends in mineral nitrogen pool size changes among models across all scenarios (Fig. 4c), indicating a large discrepancy and very limited agreement about mineral nitrogen estimation among ESMs. Better understanding and proper representation of organic nitrogen decomposition (nitrogen mineralization) are key to narrowing the gap (Thomas et al., 2015; Forsmark et al., 2020; Davies-Barnard et al., 2020).

Because of the lack of observational constraints and process-level understanding, the terrestrial carbon cycle is a major source of uncertainty contributing to future climate projections and past climate simulations (Friedlingstein et al., 2014; Ciais et al., 2014). By analyzing outputs from 12 CMIP5 ESMs, a study has found that the projection uncertainty in global land carbon storage by 2100 is >160 GtC (more than 50 % larger than that of the studied CMIP6 ESMs here), primarily driven by model structure differences (Lovenduski and Bonan, 2017). Though ESMs have been significantly improved over the past years (Eyring et al., 2019, 2021; Chen et al., 2021; Washington et al., 2009), with the inclusion of a nitrogen cycle as a realistic constraint on the carbon cycle being the most recent improvement (at the time of CMIP6) (Davies-Barnard et al., 2020; Wei et al., 2022), the spread of the carbon–concentration feedback is not much narrowed from CMIP5 to CMIP6 (Arora et al., 2020). The variability in nitrogen pool size dynamics (Fig. 4b and c) indicates considerable nitrogen-related uncertainties accompanied by the carbon–nitrogen coupling (Du et al., 2018; Thomas et al., 2013, 2015). So far, there has been no standardized nitrogen cycle validation for ESMs, largely because of limited observations, diverging representations, and the diverse upscaling approaches used in CMIP6 ESMs (Zaehle et al., 2014; Zaehle and Dalmonech, 2011; Spafford and MacDougall, 2021; Zhu et al., 2018). An improved quantification of nitrogen effects on the carbon cycle and climate, e.g., isolating the nitrogen cycle feedback from the carbon cycle feedback, though requiring meticulous experiment design and extra model simulation, might be necessary to improve climate projection uncertainty attribution (Spafford and MacDougall, 2021). The lack of nitrogen observations and limited mechanistic understanding remain fundamental challenges in reducing uncertainty related to the nitrogen effect (Zaehle et al., 2014). The forthcoming online calibration of the updated MAGICC to CMIP6 ESMs and its application (e.g., sensitivity analyses, perturbed parameter analyses, and feedback analysis) should provide new insights into the projection uncertainty, which would be too computationally expensive to obtain from ESMs (Bonan et al., 2019; Deser et al., 2020).

4.2 The nitrogen effect on NPP

The NPP simulated by CMIP6 ESMs has a consistent trend and is much more constrained than the pool sizes (Figs. 5a and A8). The NPP simulation in CMIP6 ESMs has been improved from the CMIP5 ESMs (Wei et al., 2022), attributed to advancements in nitrogen processes and the availability of more observational data (Collier et al., 2018; Randerson et al., 2009). Based on our calibrations, there is generally nitrogen limitation on global-mean, annual-mean NPP (Fig. 5b), which is consistent with experimental studies and other model simulations (LeBauer and Treseder, 2008; Wieder et al., 2015b; Thornton et al., 2007, 2009). According to our calibrations, NorESM2-LM is the only model that shows nitrogen fertilization effects on NPP during the historical period, alongside a persistent intensification of nitrogen limitation during the scenario period. This historical fertilization coincides with an enrichment of the model's mineral nitrogen pool, while the increasing limitation correlates with the accumulation of organic nitrogen (see Figs. 4b and c and A8d). These findings suggest that nitrogen mineralization is constrained by nitrogen availability during the scenario period in NorESM2-LM, resulting in reduced mineral nitrogen levels and thereby exacerbating NPP limitation.

Figure 5 (a) Net primary production (NPP) from CMIP6 ESMs and (b) the emulated nitrogen effect on NPP (ϵ_CN(NPP)) across different scenarios (the dashed blue line serves as a reference for ϵ_CN(NPP)=1).

Based on our calibrations, MPI-ESM1-2-LR and UKESM1-0-LL show the strongest nitrogen limitation on NPP during their hist_SSP runs, with an ϵ_CN(NPP) range of 0.73–0.88 and 0.77–0.90, respectively. The strong nitrogen limitation from MPI-ESM1-2-LR in its 1pctCO2 simulation matches the continuous depletion of its mineral nitrogen pool (Fig. 4c). It is noted that JSBACH, the land component of MPI-ESM1-2-LR (Mauritsen et al., 2019), shows very limited nitrogen limitation on NPP in its CMIP5 1pctCO2 idealized simulation (Goll et al., 2017). The more severe depletion of mineral nitrogen in its CMIP6 output (maximum > 0.43 GtN, Fig. 4c) than its CMIP5 result (maximum < 0.3 GtN, value from the reference publication), along with the much higher NPP simulated in CMIP6 (maximum ∼ 110 GtC yr⁻¹, Fig. 5a) compared to CMIP5 (maximum ∼40 GtC yr⁻¹, value from the reference publication), might be the reason for the strong nitrogen limitation. On the other hand, we suspect that the strong nitrogen limitation on NPP inferred for UKESM1-0-LL is primarily the result of the incongruent high-NPP and low-PU results from UKESM1-0-LL outputs (Fig. A6; i.e., this could be a calibration issue rather than an easily explained feature of UKESM1-0-LL).

Based on our calibration, CMCC-CM2-SR5 and CMCC-ESM2 exhibit the least pronounced nitrogen limitation on NPP (or even indicate fertilization) during their hist_SSP runs, with an ϵ_CN(NPP) range of 0.93–1.05 and 0.93–1.03, respectively. The similarly weak limitation is also observed in our calibration to their 1pctCO2 runs. Their organic nitrogen pool sizes keep decreasing in the 1pctCO2 experiment (Figs. 4b and A8b), which should contribute a large flux of mineral nitrogen via decomposition. However, their mineral nitrogen pools are not enriched (Figs. A6 and A8d), indicating a considerable mineral nitrogen loss from these two models. The slight nitrogen fertilization applied to NPP (ϵ_CN(NPP)>1) is found in our calibration to the CMCC models from 1975 to 2025 or 2100 depending on the scenario. The Community Land Model (CLM) serves as the land component for both CMCC models (Lovato et al., 2022). In the Duke University and Oak Ridge National Laboratory (ORNL) Free-Air CO₂ Enrichment (FACE) experiments, CLM has exhibited a nearly negligible initial-year nitrogen-based NPP response (defined as NPP $/$ canopy nitrogen), alongside a 6 %–10 % NPP response to elevated CO₂ (Zaehle et al., 2014). The near-zero nitrogen-based NPP response and significant NPP increase suggest that CLM perceives a relatively high canopy nitrogen content (as supported by the large land organic nitrogen pool sizes in the CMCC models' output in Fig. A8b). This potentially contributes to the weak nitrogen limitation observed in the CMCC models at the global scale and underscores the potential utility of regional observation studies in constraining the global nitrogen effect. Most nitrogen limitation factors fall within the range of 0.85–1.0, with limitations increasing in higher-emission scenarios. This finding aligns with the OCN results from the RCP8.5 and RCP2.6 simulations (Fig. 2 and Meyerholt et al., 2020).

Although both land surface models and our calibrations to ESMs indicate a continuous nitrogen limitation on NPP at the global scale, there is room for debate regarding the realism of a long-term nitrogen limitation. Considering the substantial amounts of atmospheric nitrogen deposition and anthropogenic additions, alongside the ubiquitous presence of nitrogen-fixing organisms, the ability of ecosystems to optimize nitrogen use efficiency in response to varying nitrogen availability remains unclear. Previous studies suggest that ecosystems should be capable of effectively balancing these factors to alleviate long-term limitations on overall NPP (Vitousek and Howarth, 1991). This holds particularly true for tropical forests where nitrogen is abundant and rapidly circulated (Hedin et al., 2005; Cusack et al., 2011). In such cases, other nutrients like phosphorus and potassium might play a critical role (Wright et al., 2011; Alvarez-Clare et al., 2013). Further validation is required to assess the long-term effects of nitrogen limitation on NPP and to understand the differences between regional and global patterns.

4.3 The nitrogen effect on pool turnovers

According to our calibrations, the carbon–nitrogen coupling has largely enhanced plant carbon turnover in both CMCC models and MIROC-ES2L, as evidenced by ${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LP}}_{\mathrm{c}}\right)}$ values ranging from 1.09 to 1.50, 1.07 to 1.50, and 1.12 to 1.40, respectively, during the hist_SSP period (Fig. 6a). The ${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LP}}_{\mathrm{c}}\right)}$ values exhibit an increasing trend from the SSP1-2.6 scenario to the SSP5-8.5 scenario. The calibration to the other ESMs, however, suggests that the plant nitrogen status has inhibited litter production (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LP}}_{\mathrm{c}}\right)}$ in the range of 0.6–0.8 during the hist_SSP period and no significant difference among different scenarios). The carbon–nitrogen coupling effect (ϵ_CN) and temperature effect (ϵ_dT) effectively change the turnover rate determined by the initial turnover times (Eq. 29). Considering that the temperature feedback is always 1 at the beginning of each experiment (dT=0) and the dynamics of plant carbon are similar across different ESMs (Figs. 4a and A6), the strong nitrogen enhancement in litter carbon production in the CMCC models is mainly because of their high initial plant carbon turnover times (31 and 66 years, respectively, Table A2), compared to the other ESMs (15 to 25 years, Table A2). The combination allows for the total litter production rate to remain at a similar level for all the ESMs, a key requirement given that they have similar plant pool size changes (Fig. A6), relatively consistent NPP (the influx for the plant carbon pool, Fig. 5), and similar plant litter respiration (the outflux for the plant carbon pool, the LPR₀ parameter in Table A2). The enhanced litter production in MIROC-ES2L, however, is needed to compensate for its much lower plant respiration flux (LPR₀=0.96 GtC yr⁻¹ in MIROC-ES2L vs. 4.16–9.90 GtC yr⁻¹ in other ESMs). The weak limitation (or even fertilization) of NPP found for the CMCC models (Fig. 5b) and the continuous loss of soil organic nitrogen (Figs. 4b and A6, 1pctCO2 and historical period) suggest that the system is less nitrogen limited. The resulting plant nitrogen availability partially contributes to the fast plant carbon turnovers in our calibration to the CMCC models.

Figure 6 The emulated nitrogen effect on carbon pool turnovers, including (a) litter production (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LP}}_{\mathrm{c}}\right)}$), (b) litter decomposition (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LD}}_{\mathrm{c}}\right)})$, and (c) soil respiration (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{SR}}_{\mathrm{c}}\right)}$), from CMIP6 ESMs across different scenarios.

Our calibrations to all the ESMs except for MPI-ESM1-2-LR show consistently inhibited litter decomposition after the nitrogen effect is applied (${\mathit{ϵ}}_{CN\left({\mathrm{LD}}_{\mathrm{c}}\right)}$ in the range of 0.55–0.96 for the hist_SSP runs), and such inhibition slightly increases from SSP1-2.6 to SSP5-8.5. The significantly enhanced litter decomposition in our calibration to MPI-ESM1-2-LR is attributed to the high carbon input into its litter pool and the small temperature sensitivity, which are supported by (1) the highest NPP and its partition into the litter pool among the ESMs (Fig. 5a and Table A2, f_NPP2L=0.59), (2) the highest litter production partition into the litter pool ($f_{\mathrm{LP}\mathrm{2}{\mathrm{L}}_{\mathrm{c}}}=\mathrm{0.96}$, Table A2), and (3) the lowest temperature sensitivity of litter decomposition ($s_{\mathrm{d}T\left({\mathrm{LD}}_{\mathrm{c}}\right)}=\mathrm{0.024}$, Table A2). The strong nitrogen limitation on NPP (Fig. 5b), the continuous depletion of the mineral nitrogen pool (Figs. 4c and A6, 1pctCO2 and historical period), and the highest biological nitrogen fixation (Fig. 7a) in our calibration to MPI-ESM1-2-LR indicate that the plant nitrogen is insufficient. The strong ${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{LD}}_{\mathrm{c}}\right)}$ along with the low $f_{\mathrm{LD}\mathrm{2}{\mathrm{S}}_{\mathrm{c}}}$ (0.07, Table A2) suggests that the system is trying to mineralize more litter carbon to mediate the plant nitrogen deficiency and to maintain an ecologically reasonable carbon : nitrogen stoichiometry.

Figure 7 (a) Biological nitrogen fixation (BNF) and (b) its ratio to net primary production (BNF : NPP) from CMIP6 ESMs across different scenarios.

The soil respiration is found to be restricted in our calibration to UKESM1-0-LL (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{SR}}_{\mathrm{c}}\right)}=\mathrm{0.61}$–0.69 during the hist_SSP period), while it is significantly enhanced in other ESMs. The strongest enhancement is found in the two CMCC models (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{SR}}_{\mathrm{c}}\right)}=\mathrm{1.24}$–1.50 during the hist_SSP period), followed by NorESM2-LM (${\mathit{ϵ}}_{\mathrm{CN}\left({\mathrm{SR}}_{\mathrm{c}}\right)}=\mathrm{1.20}$–1.36 during the hist_SSP period). The calibrations to these three models show longer soil carbon turnover times (283–476 years, Table A2) than the other ESMs. The order-of-magnitude-larger mineral nitrogen pool size in NorESM2-LM compared to other ESMs (Fig. A8d) and the continuously growing mineral nitrogen and organic nitrogen pool sizes (Fig. 4c, hist_SSP experiments) support the enhanced soil organic matter decomposition. It is observed that the mineral nitrogen pool size exhibits a continuous decrease during the NorESM2-LM's 1pctCO2 run, while it first increases and then decreases in its SSP5-8.5 run (Figs. 4d and A6). Considering that plant uptake and net mineralization are the two major fluxes controlling the mineral nitrogen dynamics, this result suggests a potential threshold associated with climate or CO₂ concentration, limiting the net mineralization rate from matching the ongoing increase in plant uptake (Fig. A6). The high temperature change (Fig. A7) and its subsequent large temperature effect on respiration in UKESM1-0-LL could be responsible for the nitrogen-inhibited soil respiration. The substantial land carbon accumulation in MIROC-ES2L (Fig. 4a) requires less respiration, thus explaining the neglectable nitrogen effect on its soil respiration. The diverse impacts of nitrogen on soil carbon turnover align with existing experimental findings, which have demonstrated contrasting trends in nitrogen additions across various substrate decompositions (Averill and Waring, 2018; Hobbie, 2008). As a result, studies have suggested that the classic stoichiometric decomposition theory should be revised (Craine et al., 2007).

5 Limitations and implications

A key limitation of RCMs is their resolution, which arises from the trade-off between spatial heterogeneity and computational efficiency. Therefore, it is important to note that the feedbacks, pool size dynamics, and fluxes discussed in this paper are aggregated from extensive spatiotemporal datasets into an annual and global framework. This means the representation of the carbon–nitrogen cycle reflects a synthesis of diverse regional and sub-annual dynamics. Consequently, the results presented here are at a global scale and may differ from regional or sub-regional studies, necessitating cautious interpretation.

5.1 The simulation of mineral nitrogen pool dynamics

Our emulation is less effective in capturing the mineral nitrogen pool sizes for some ESMs (e.g., NorESM2-LM, Figs. 3 and A6). However, the large mineral nitrogen pool in NorESM2-LM (Fig. A8d) and the significant discrepancies in mineral nitrogen pool size changes among ESMs (Fig. 4c) and land surface models (Fig. 2) highlight the need for a deeper understanding of these dynamics in complex models. This should include enhanced theoretical underpinning and improved observational constraints (Zaehle and Dalmonech, 2011; Thomas et al., 2015). The heterogeneous terrestrial nitrogen cycle results in challenges and discrepancies even compared to the regional observations we do have (Schulte-Uebbing and de Vries, 2018; Ramm et al., 2022; Menge et al., 2012). The uncertain atmospheric deposition and nitrogen fertilizer application further complicate the evaluation of the mineral nitrogen pool size and its dynamics (Mulvaney et al., 2009; Reay et al., 2008; Gruber and Galloway, 2008).

The two main controls of the mineral nitrogen pool size, the nitrogen mineralization (from organic decomposition) and mineral nitrogen loss, are still poorly understood at the process level (Manzoni et al., 2008; Hedin et al., 2005). For instance, the microbial decomposition of organic matter (heterotrophic respiration) can be limited, stimulated, or even unaffected by the nitrogen addition as a result of differences in soil microbial biomass or activity changes (Bardgett et al., 1999). The nitrogen effect on decomposition has been found to be sensitive to the types of substrates, but generally the impact on the decomposition rate is negative or neutral (Hobbie, 2008). The root allocation, plant growth, litter production, biodiversity, etc., are all influenced by nitrogen (Phoenix et al., 2006; Wright et al., 2011). However, a 13-year-long nitrogen addition study has found that lower nitrogen addition rates had no effect on litter production or soil respiration in a Pinus sylvestris forest (Forsmark et al., 2020). This raises questions about the overall nitrogen effect at the ecosystem level, particularly considering the uneven geographical distribution of atmospheric nitrogen deposition (Phoenix et al., 2006). Alongside climate factors such as warming and precipitation, as well as other ecological or physical constraints, the situation becomes even more complicated (Plett et al., 2020; Lim et al., 2015; Cusack et al., 2010; Reich et al., 2014; Li et al., 2019). These findings underscore the highly complex dynamics of mineral nitrogen, suggesting that the current formulation of mineral nitrogen loss in CNit, characterized by simple first-order decay and a single temperature response (Eq. 29), is very likely an oversimplification. Introducing constraints into the global stoichiometry of nitrogen mineralization could potentially enhance the modeling of mineral nitrogen pool dynamics (Manzoni et al., 2008; Meyerholt and Zaehle, 2015). Nevertheless, these new developments in MAGICC are a step towards better representation of carbon–nitrogen dynamics. The calibration results also demonstrate that the current CNit formulation mostly captures the trend of mineral nitrogen pool sizes in the majority of the studied CMIP6 ESMs and experiments at the global-mean, annual-mean level (Figs. 3 and A6).

5.2 The biological nitrogen fixation as an input instead of being simulated

Biological nitrogen fixation (BNF) serves as the primary non-anthropogenic nitrogen input in the global nitrogen cycle (Vitousek et al., 2002; Gruber and Galloway, 2008; Fowler et al., 2013). The trend of BNF generally mirrors that of NPP (Figs. 7a and 5a), but the relative differences between ESMs are more pronounced (e.g., the largest initial biological nitrogen fixation found in MPI-ESM1-2-LR is approximately 4 times that of the smallest one found in NorESM2-LM). This similarity in the trend arises because the model representations of BNF in all the studied ESMs were based, at least partly, on its empirical relationship with either NPP or evapotranspiration derived from the widely recognized global BNF estimate (Cleveland et al., 1999). Figure 7b illustrates a variety of global patterns of BNF : NPP in the CMIP6 ESM outputs, showcasing diverse trends such as decreases (CMCC-CM2-SR5, CMCC-ESM2, and MIROC-ES2L in the 1pctCO2 experiment), increases (UKESM1-0-LL in the 1pctCO2 experiment and MIROC-ES2L in SSP experiments), stability (UKESM1-0-LL in SSP scenarios), or even a peak followed by a decrease (NorESM2-LM in low-SSP scenarios and MIROC-ES2L in high-SSP scenarios). It is noteworthy that the mean BNF : NPP remains relatively constant at approximately 0.0018 GtN GtC⁻¹ across all experiments.

Recent findings from a meta-analysis of field measurements challenge the notion of a statistically significant relationship between biological nitrogen fixation and NPP/evapotranspiration (Davies-Barnard and Friedlingstein, 2020). An analysis of model uncertainty in recent studies highlighted that variations in the BNF representation could significantly impact future climate projections (Wieder et al., 2015a; Kou-Giesbrecht and Arora, 2022). For instance, employing different BNF representations within a shared framework has resulted in modeled BNF responses to elevated (200 ppm higher) CO₂ ranging from $-\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{56}\times {\mathrm{10}}^{-\mathrm{3}}$ GtN yr⁻¹ (Meyerholt et al., 2016). This variation has led to a global land carbon storage range of 281 to 353 GtC (over ∼150 years of simulation), with N₂O emissions fluctuating from $-\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{3}}$ GtN yr⁻¹ (Meyerholt et al., 2016). A recent study assessing BNF structural uncertainty in CMIP6 ESMs has revealed that the response of BNF and other nitrogen cycle variables could differ, even among models with similar structures (Davies-Barnard et al., 2022). Conflicting empirical relationships and updated observations underscore the considerable uncertainty in BNF estimation and potential need for revisions to BNF formulations. While we could develop parameterizations to emulate the BNF formulations used in ESMs, the BNF flux represents a relatively minor flux with questionable data quality and highly uncertain formulations and/or mechanisms (e.g., CABLE prescribes constant BNF). Therefore, we opt not to pursue it further here. Instead, CNit directly prescribes BNF from CMIP6 ESM outputs to circumvent further structural uncertainty stemming from simplified parameterization.

5.3 The disentangled climate feedback and nitrogen effect from emulation

One advantage of RCMs is that their simplified formulations attempt to capture the overall effects of complex processes in ESMs, aiding the identification and quantification of key effects in the system. Based on the assumptions and definitions in this updated carbon–nitrogen cycle, we can separate the temperature feedback and carbon–nitrogen coupling feedback for different pool turnovers (Eq. 29). However, a significant challenge in this separation arises from the exponential formulation of both feedbacks and the increasing trends of all the feedback proxies (temperature change (dT), nitrogen plant uptake (PU), and nitrogen atmospheric deposition (AD), Figs. A4 and A6). This setup suggests that these feedbacks (or the related sensitivity parameters) may change in opposite directions to compensate each other while still producing a similar overall feedback. The formulations (Eq. 29) suggest that the turnover time and the overall feedback can also offset each other to reach a similar turnover flux.

To examine the correlation of parameter values and feedback separation, we applied Markov chain Monte Carlo (MCMC) sampling to the sensitivity parameters and turnover times for each of the individual ESMs (60 walkers × 1000 iterations = 60 000 runs, starting from the best-estimate parameter values). The results show that there is a weak-to-moderate negative correlation between temperature sensitivity and plant uptake sensitivity for most of the ESMs (the absolute value of Spearman's r=0.04–0.72 for most cases, highlighted in yellow in Fig. A9). The strongest correlations are found in NorESM2-LM (Spearman's $r=-\mathrm{0.53}$, −0.72, and −0.46 for the litter production, litter decomposition, and soil respiration, respectively). The weakest correlations are found in MIROC-ES2L (Spearman's $r=-\mathrm{0.04}$, −0.34, and 0.04 for the litter production, litter decomposition, and soil respiration, respectively). The temperature sensitivity and atmospheric deposition sensitivity show relatively weak correlations (the absolute value of Spearman's r<0.3 for most cases, Fig. A9).

The weak-to-moderate correlations between temperature sensitivity and plant uptake sensitivity are mainly due to the feedback proxies, dT and PU, which, though both exhibiting the same increasing trend, do not strictly change with the same gradients (e.g., the temperature change fluctuates near zero, while plant uptake shows a clear increasing trend from 1850–1975 in the historical simulations, Figs. A4 and A6). Increasing one type of sensitivity while decreasing another, though it could lead to similar overall feedback at the early stage when temperature change and plant uptake are less perturbed, cannot guarantee similar overall feedback throughout the entire time series (e.g., the difference in overall feedback from different sensitivities is amplified as the temperature change gets higher and plant uptake becomes larger). In other words, the parameter values cannot simply all vary in the opposite direction to compensate for the feedback. Instead, to reach the desirable turnover flux to satisfy the pool size dynamics, the parameter values need to be adjusted (and not necessarily offset each other) to obtain the respective “correct” or best-estimate feedback. Such results suggest that separating the carbon–nitrogen coupling feedback from the temperature feedback (the assumption for Eq. 29) is a reasonable assumption, although this should be investigated further in future work and is achieved most readily in simulations that specifically target this separation, for example the nitrogen-on and nitrogen-off experiments performed by the land surface models we used here.

The turnover times and plant uptake sensitivities exhibit strong positive correlations in all the ESMs (displayed in bold in Fig. A9), indicating that plant uptake sensitivity is the predominant factor influencing the overall feedback to compensate for turnover time changes. This is supported by the dominance of carbon–nitrogen coupling feedback in both the magnitudes and the trends of the overall feedback (Figs. 6 and A10), further emphasizing the substantial disparity between temperature feedback and carbon–nitrogen coupling feedback and the imperative of distinguishing between the two.

One limitation of disentangling the carbon–nitrogen coupling feedback from the climate feedback is that the feedback strength is primarily derived from emulation. Although several factors support this distinction, (1) the evidence presented in Figs. A9 and A10 underscores the clear differentiation between climate feedback and carbon–nitrogen coupling feedback, (2) the dynamics of the pool size offer indirect yet compelling constraints (referred to as “emergent” constraints) for the feedback, and (3) the selection of feedback formulations (exponential rather than linear relationships) and proxies (with varying magnitudes) further restricts the parameters from offsetting each other. However, the absence of nitrogen-off simulations from the CMIP6 ESMs presents challenges for direct verification. Given the computational expense of running all scenarios in nitrogen-off mode, it is recommended that ESMs perform nitrogen-off simulations for select idealized scenarios (e.g., 1pctCO2 or flat10 – constant emissions of CO₂ of 10 GtC yr⁻¹; Sanderson et al., 2024) for diagnostic purposes.

6 Conclusion and future work

In this work, we have detailed CNit v1.0, a new coupled carbon–nitrogen model for MAGICC, which enhances the model's capability to represent terrestrial biogeochemical processes and feedbacks. Based on the offline calibration results from land surface models and multiple CMIP6 ESMs, we have demonstrated that CNit v1.0 is able to effectively emulate the behavior of the carbon–nitrogen cycles from various more complex models, encompassing a broad spectrum of carbon–nitrogen states and dynamics.

The temperature change and carbon–nitrogen state/dynamics (especially those related to the nitrogen cycle and mineral nitrogen) exhibit significant variability among CMIP6 ESMs, particularly in their 1pctCO2 and high-SSP scenario runs, which highlights the model structure uncertainty. The contrasting trends in mineral nitrogen dynamics and the order-of-magnitude differences in pool sizes underscore the limited agreement in mineral nitrogen modeling. A thorough analysis, focusing on the new uncertainties introduced by the nitrogen cycle, is imperative for the CMIP6 and future ESMs. The upcoming (in future research) sensitivity analysis, perturbation parameter analysis, and feedback analysis of the updated MAGICC model are expected to provide insights for uncertainty attribution.

The CNit emulation indicates a general nitrogen limitation on NPP, which follows a similar trend across the studied CMIP6 ESMs. Combining the results from NPP and turnovers suggests that, at the multi-model mean level, the carbon–nitrogen coupling limits both NPP and plant and litter carbon pool turnovers, though the weaker NPP nitrogen limitation could also lead to significantly enhanced litter production. Soil respiration is instead enhanced in most of the ESMs. The combination indicates that terrestrial ecosystems may become net carbon sources sooner than we would expect based on models that do not consider the impact of the nitrogen cycle.

The presented carbon–nitrogen coupling in MAGICC demonstrates the ability to emulate many complex models while nonetheless having limitations, particularly in simulating mineral nitrogen pool dynamics and biological nitrogen fixation. There are currently significant inconsistencies between ESM outputs and observations of the mineral nitrogen pool size and biological nitrogen fixation, in terms of both magnitudes and trends, suggesting that substantial revisions are possible in the near future. Therefore, the current formulation and treatment of these aspects in MAGICC may have to be updated too while aiming to continue to strike a balance between model simplicity, process representation, and emulation performance, reflecting a fundamental design principle for RCMs and MAGICC in particular. Future work on MAGICC's carbon–nitrogen cycle will focus on the calibration of the full MAGICC structure to CMIP6 ESMs (and/or observational data), evaluation of model performance with respect to computational efficiency and mechanistic insight, incorporation of additional constraints, uncertainty quantification, sensitivity analysis, application of probabilistic projections, and continued model development (e.g., land use emission implementation and nitrogen process representation) to align with advances in complex models and emerging theoretical frameworks.

Appendix A

Table A1 Full list of calibrated CNit parameters: long name and range.

Figure A1 Illustration of the functionality of the method factor for CO₂ fertilization ($m_{{\mathrm{CO}}_{\mathrm{2}}}$) and temperature feedback (m_dT). An $m_{{\mathrm{CO}}_{\mathrm{2}}}$ of 0, 1, and 2 represents the logarithmic, rectangular, and sigmoidal CO₂ fertilization formulations, respectively (Eq. 17). Similarly, an m_dT of 0 and 1 corresponds to the exponential and sigmoidal temperature response formulations, respectively (Eq. 20). Intermediate values represent a linear combination of the two formulations.

Figure A2 Relationship between nitrogen plant uptake and net mineralization as simulated by CABLE and OCN.

Figure A3 Global average surface temperature change over land (dT_LAND, delta annual mean tas over land) and CO₂ concentrations from land surface models across different scenarios.

Figure A4 Global average surface temperature change over land (dT_LAND, delta annual mean tas over land) and CO₂ concentrations from CMIP6 ESMs across different scenarios.

Table A2 Full list of calibrated CNit parameters: values.

Note that the missing values (indicated by “–”) for UKESM1-0-LL are due to the model's lack of a litter pool. As a result, there are no turnover times, feedback-related parameters, or fractionation values associated with the litter pool.

Figure A5 Comparison of net primary production (NPP), plant carbon pool size (C_P), litter carbon pool size (C_L), soil carbon pool size (C_S), land carbon pool size (C_LAND), net land carbon flux (fC_NetLAND), nitrogen plant uptake (PU), plant nitrogen pool size (N_P), litter nitrogen pool size (N_L), soil nitrogen pool size (N_S), mineral nitrogen pool size (N_M), land organic nitrogen pool size (N_LANDo, sum of nitrogen in plant, litter, and soil pools), and net land nitrogen flux (fN_NetLAND) between CABLE or OCN outputs (blue lines) and CNit emulations (orange lines). The experiments labeled C, CN, and CNd denote the carbon-only, carbon–nitrogen coupled with constant nitrogen atmospheric deposition, and carbon–nitrogen coupled with dynamic nitrogen atmospheric deposition configurations, respectively, in the land surface models.

Figure A6 Comparison of net primary production (NPP), plant carbon pool size (C_P), litter carbon pool size (C_L), soil carbon pool size (C_S), land carbon pool size (C_LAND), net land carbon flux (fC_NetLAND), nitrogen plant uptake (PU), plant nitrogen pool size (N_P), litter nitrogen pool size (N_L), soil nitrogen pool size (N_S), mineral nitrogen pool size (N_M), land organic nitrogen pool size (N_LANDo, sum of nitrogen in plant, litter, and soil pools), and net land nitrogen flux (fN_NetLAND) between CMIP6 ESM outputs (blue lines) and CNit emulations (orange lines).

A1 The diversity of temperature outputs from CMIP6 ESMs

Temperature change, a pivotal driving force for the carbon–nitrogen cycle, exhibits significant variation between the two land surface models (Fig. A3) and among the simulations of CMIP6 ESMs (Fig. A4), even when they undergo the same experiment. For the sake of comparison, Fig. A7 and the subsequent discussions focus on the common experimental periods for scenarios (e.g., 1850–2100 for the hist_SSP experiments).

UKESM1-0-LL shows the highest temperature change among all models and experiments, whereas NorESM2-LM exhibits the lowest temperature change. Both the idealized 1pctCO2 – one of the base experiments in the Diagnostic, Evaluation and Characterization of Klima (DECK) experiments – and the historical simulation are in the core set of experiments performed under CMIP5, CMIP6, and previous CMIP phases (Eyring et al., 2016; Taylor et al., 2012). As CMIP6 and C4MIP necessitate consistent forcings and experimental protocols for simulations conducted by participating ESMs (Eyring et al., 2016; Jones et al., 2016), the wide spread of the land surface temperature change values – especially from the 1pctCO2 experiment where land use change is not included (with a standard deviation of 1.3 °C and absolute difference of 3.9 °C at the end of the simulation) – highlights the various parameterizations of physical processes in ESMs resulting in large differences in the ESMs climate sensitivities (Rugenstein et al., 2020; Meehl et al., 2020), for example, the structural uncertainty (Deser et al., 2020; Duan et al., 2021). The previous MAGICC simulation with constraints from historical CO₂ measurements and temperature observations is found to reduce uncertainty in the temperature projections (Bodman et al., 2013).

Recent studies interpreting surface air temperature outputs from multiple CMIP6 ESMs indicate that the multi-model mean effectively captures the historical temperature trend in observations (Fan et al., 2020; Papalexiou et al., 2020). Results from a study using outputs from 29 CMIP6 ESMs show that the post-1988 warming is overestimated in 90 % of the simulations and the observed long-term persistence of global mean temperature (for the period of 1880–2014) is not accurately captured in most of the ESMs (Papalexiou et al., 2020), suggesting further model selections based on the case-specific intended uses. However, previous evaluation of the long-term persistence of temperature on continental areas (60° S–60° N) during 1930–2004 in CMIP5 ESMs demonstrated that most models captured the long-term persistence reasonably well (Kumar et al., 2013). Moreover, grouping CMIP6 ESMs and re-analyzing the global mean temperature based on the grouped models can also lead to different conclusions on the warming trend (Scafetta, 2023). These results indicate that more careful interpretation of the simulated temperature is needed. They also justify using each ESM's global mean land temperature as input in this study instead of the global mean temperature (to avoid differences in calibration based on inconsistency with the target model's temperature rather than any issue with the reduced-complexity model).

Figure A7 Global average surface temperature change over land (dT_LAND, delta annual mean tas over land) from CMIP6 ESMs across different scenarios (over the common experimental period).

A2 The diversity of the carbon–nitrogen cycle in CMIP6 ESMs

Based on the varied temperature results, it is not surprising that the carbon–nitrogen cycle fluxes and pools from the CMIP6 ESMs are diverse (Fig. A8). The initial land carbon pool ranges from 1396 GtC (MPI-ESM1-2-LR) to notably higher values of 3300 GtC (CMCC-CM2-SR5 and CMCC-ESM2). The initial nitrogen pool sizes are even more inconsistent among ESMs. The largest initial organic nitrogen pool is 265 GtN from the two CMCC models, which is more than 4 times the size of the smallest one from MPI-ESM1-2-LR (61 GtN). These results have led to a wide range of initial organic carbon : nitrogen ratios from 12 (the two CMCC models) to 23 (MPI-ESM1-2-LR) (Fig. A8c). The trends for the carbon pool size and carbon : nitrogen ratio exhibit a similar pattern. They display a consistent increase in the 1pctCO2 scenario, while in the hist_SSP simulations, they initially decrease (1850–1970) and then rise again (1970–2100). The mineral nitrogen pool, on the contrary, shows significant variations in both pool sizes and trends across CMIP6 ESMs (Fig. A8d). The initial pool sizes range from <0.2 GtN (UKESM1-0-LL and MIROC-ES2L) to 5.6 GtN (NorESM2-LM). The trends are found to be either opposite or unrelated.

Figure A8 Diversity of land carbon pool size (C_LAND), land organic nitrogen pool size (N_LANDo, sum of nitrogen in plant, litter, and soil pools), carbon : nitrogen ratio (CN ratio), and mineral nitrogen pool size (N_M) from CMIP6 ESM outputs.

The substantial variation in simulated carbon pools is a long-standing issue for both CMIP5 and CMIP6 ESMs (Anav et al., 2013; Varney et al., 2022). The initial-condition differences are responsible for the models' internal variability (Deser et al., 2020; Kumar and Ganguly, 2018), which accounts for more than half of the inter-model spread in near-term climate projections (Deser et al., 2012). Such differences also contribute to their respective carbon cycle projections. The different initial carbon–nitrogen cycle state among different ESMs further complicates their comparison (Spafford and MacDougall, 2021). The standard deviation of initial land carbon pool sizes is different for the 1pctCO2 (699 GtC) and historical (720 GtC) scenarios. Four of the studied ESMs have used nearly the same starting land carbon pool size for both scenarios (difference < 1.5 GtC), while UKESM1-0-LL and NorESM2-LM have a large difference of 9 and 184 GtC (out of a total of roughly 2300 GtC in UKESM1-0-LL and 3000 GtC in NorESM2-LM), respectively. The varying initial pool sizes can pose significant challenges for emulators employing first-order decay for pool turnovers as the turnover time predominantly influences the magnitude of the “base” turnover flux (e.g., without any feedback scalers). Since MAGICC emulation has used the same set of parameters (including the turnover times) to emulate all the experiments, it explains the jump in the emulated soil carbon pool sizes in the 1pctCO2 experiment for these two models (Fig. A6).

Figure A9 Correlation of turnover times and feedback-related parameters from the CMIP6 ESMs. The numbers indicate Spearman's correlation coefficients (r) between pairs of parameters, with ^* denoting p values < 0.001. Correlations between temperature sensitivities and plant nitrogen uptake sensitivities are highlighted in yellow, while correlations between turnover times and plant nitrogen uptake sensitivities are shown in bold. Missing values for UKESM1-0-LL are due to the absence of a litter pool in this model, resulting in no turnover time or feedback-related parameters for the litter pool.

Figure A10 The emulated overall and temperature effect on carbon pool turnovers, including litter production (${\mathit{ϵ}}_{\mathrm{d}T\mathit{\&}\mathrm{CN}\left({\mathrm{LP}}_{\mathrm{c}}\right)}$ and ${\mathit{ϵ}}_{\mathrm{d}T\left({\mathrm{LP}}_{\mathrm{c}}\right)}$, litter decomposition (${\mathit{ϵ}}_{\mathrm{d}T\mathit{\&}\mathrm{CN}\left({\mathrm{LD}}_{\mathrm{c}}\right)}$ and ${\mathit{ϵ}}_{\mathrm{d}T\left({\mathrm{LD}}_{\mathrm{c}}\right)}$), and soil respiration (${\mathit{ϵ}}_{\mathrm{d}T\mathit{\&}\mathrm{CN}\left({\mathrm{SR}}_{\mathrm{c}}\right)}$ and ${\mathit{ϵ}}_{\mathrm{d}T\left({\mathrm{SR}}_{\mathrm{c}}\right)})$ from CMIP6 ESMs across different scenarios.