Evaluating nitrogen cycling in terrestrial biosphere models: a disconnect between the carbon and nitrogen cycles

Terrestrial carbon (C) sequestration is limited by nitrogen (N), an empirically established constraint that could intensify under CO₂ fertilization and future global change. The terrestrial C sink is estimated to currently sequester approximately a third of annual anthropogenic CO₂ emissions based on an ensemble of terrestrial biosphere models, which have been evaluated in their ability to reproduce observations of the C, water, and energy cycles. However, their ability to reproduce observations of N cycling and thus the regulation of terrestrial C sequestration by N have been largely unexplored. Here, we evaluate an ensemble of terrestrial biosphere models with coupled C–N cycling and their performance at simulating N cycling, outlining a framework for evaluating N cycling that can be applied across terrestrial biosphere models. We find that models exhibit significant variability across N pools and fluxes, simulating different magnitudes and trends over the historical period, despite their ability to generally reproduce the historical terrestrial C sink. Furthermore, there are no significant correlations between model performance in simulating N cycling and model performance in simulating C cycling, nor are there significant differences in model performance between models with different representations of fundamental N cycling processes. This suggests that the underlying N processes that regulate terrestrial C sequestration operate differently across models and appear to be disconnected from C cycling. Models tend to overestimate tropical biological N fixation, vegetation C : N ratio, and soil C : N ratio but underestimate temperate biological N fixation relative to observations. However, there is significant uncertainty associated with measurements of N cycling processes given their scarcity (especially relative to those of C cycling processes) and their high spatiotemporal variability. Overall, our results suggest that terrestrial biosphere models that represent coupled C–N cycling could be overestimating C storage per unit N, which could lead to biases in projections of the future terrestrial C sink under CO₂ fertilization and future global change (let alone those without a representation of N cycling). More extensive observations of N cycling processes and comparisons against experimental manipulations are crucial to evaluate N cycling and its impact on C cycling and guide its development in terrestrial biosphere models.

1 Introduction

The terrestrial biosphere is estimated by the Global Carbon Project (GCP) to currently sequester approximately a third of anthropogenic CO₂ emissions (Friedlingstein et al., 2022). The GCP annually reports an estimate of the global carbon (C) budget, which includes an estimate of the atmosphere–land CO₂ flux based on simulations of an ensemble of terrestrial biosphere models – the trends in the land carbon cycle project (TRENDY) ensemble. In recent years, the majority of the models within the TRENDY ensemble have incorporated a representation of coupled C and nitrogen (N) cycling given the empirically established importance of N limitation of vegetation growth (Elser et al., 2007; LeBauer and Treseder, 2008; Wright et al., 2018). Whereas only 4 out of 9 models represented coupled C–N cycling in the 2013 GCP, 11 out of 16 models represented coupled C–N cycling in the 2022 GCP (Fig. 1). Capturing N constraints on C cycling is critical for realistically simulating the terrestrial C sink, which arises from the combined effects of concurrently acting global change drivers that are each modulated by N. CO₂ fertilization is limited by N (Terrer et al., 2019; S. Wang et al., 2020), intensifying N deposition increases N supply (O'Sullivan et al., 2019; Wang et al., 2017), rising temperature and varying precipitation modulate decomposition and soil N availability (Liu et al., 2017), and land use change and associated N fertilization regimes determine N supply to crops.

Figure 1 Number of terrestrial biosphere models contributing to the Global Carbon Project (the TRENDY ensemble) with and without coupled C–N cycling.

The TRENDY ensemble has been extensively evaluated against observations of the C, water, and energy cycles (Collier et al., 2018; Friedlingstein et al., 2022; Seiler et al., 2022). Within the GCP itself, the primary simulated C pools, C fluxes, and water fluxes are evaluated using a skill score system developed by the International Land Model Benchmarking (ILAMB) project that quantifies model performance by comparing model simulations to observations (Collier et al., 2018; Friedlingstein et al., 2022). ILAMB scores encompass the mean and variability of a given variable (pool or flux) over monthly to decadal temporal scales and over grid cell to global spatial scales. However, despite its importance in regulating C cycling, N cycling has not been explicitly evaluated. This is in part due to the relatively recent incorporation of N cycling in terrestrial biosphere models (Fig. 1) (Fisher and Koven, 2020; Hungate et al., 2003) but also due to the paucity of global observation-based datasets of N cycling. N exists in many forms and is lost from terrestrial ecosystems via numerous pathways (emissions of NH₃, N₂O, NO_x, and N₂, as well as NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ leaching), N processes are generally not measured in situ in networks such as FLUXNET, and remote sensing methodologies for measuring N processes are still in their infancy. Additionally, N processes exhibit extremely high spatial and temporal variabilities and are thus challenging to measure. As such, N cycling has commonly been evaluated by comparing simulated N pools and fluxes to global totals based on a small number of observations that have been scaled up or averaged to yield a value with wide confidence intervals (Davies-Barnard et al., 2020).

N cycling is implicitly evaluated by comparing terrestrial biosphere models without N cycling to those with coupled C–N cycling in reproducing observations of the C, water, and energy cycles in the absence of N cycle observations. Results suggest that there are only minor differences between the performance of models with and without N cycling. There is no significant difference between the terrestrial C sink simulated by the TRENDY models with and without N cycling (Friedlingstein et al., 2022) or between the terrestrial C sink simulated by the models participating in the Multi-scale synthesis and Terrestrial Model Intercomparison Project (MsTMIP) with and without N cycling (Huntzinger et al., 2017). Comparing the mean score across all C, water, and energy cycle variables between TRENDY models with and without N cycling yielded no significant difference (Seiler et al., 2022). However, TRENDY models without N cycling had significantly higher scores for net biome productivity than TRENDY models with N cycling (although all other variables were not significantly different between TRENDY models with and without N cycling, including vegetation C, soil C, net biome productivity, leaf area index, latent heat flux, and runoff) (Seiler et al., 2022). Despite this seeming absence of a difference between models with and without coupled C–N cycling in simulating the current terrestrial C sink, it is imperative that N constraints on C cycling are properly represented by terrestrial biosphere models in order to realistically simulate the terrestrial C sink under future global change, which modifies the C–N balance through N limitation of CO₂ fertilization and intensifying N deposition among other effects of global change. As such, explicitly evaluating N cycling processes themselves is necessary to assess the ability of terrestrial biosphere models to capture the underlying mechanisms that determine terrestrial C sequestration and thus to realistically project the future terrestrial C sink under global change.

Here, we synthesize the N pools and fluxes simulated by 11 terrestrial biosphere models in the TRENDY ensemble that participated in the 2022 GCP. We evaluate their performance in reproducing observations of three key variables of the N cycle: biological N fixation, vegetation C : N ratio, and soil C : N ratio. These three variables are critical to C cycling because (1) biological N fixation is the dominant natural N supply to terrestrial ecosystems, influencing the degree of N limitation of plant growth and thus terrestrial C sequestration, and (2) vegetation and soil C : N ratios reflect assimilated C per unit N and thus terrestrial C storage.

2 Methods

2.1 Simulation protocol

For the 2022 GCP (version 11), the TRENDY ensemble consisted of 16 terrestrial biosphere models, 11 of which represent N cycling (CABLE-POP, CLM5.0, DLEM, ISAM, JSBACH, JULES-ES, LPJ-GUESS, LPX-Bern, OCNv2, ORCHIDEEv3, and SDGVM). Although SDGVM includes a representation of N cycling, its representation is simplistic and was therefore not included. Additionally, CLASSIC contributed to the 2022 GCP without coupled C–N cycling; the S3 simulation from the TRENDY protocol was repeated by CLASSIC with coupled C–N cycling following the 2022 GCP protocol and was used here. Overall, we analysed 11 models with coupled C–N cycling (Table 1).

Table 1 Terrestrial biosphere models in the TRENDY-N ensemble and descriptions of their representations of N limitation of vegetation growth, biological N fixation, vegetation response to N limitation (i.e., strategies in which vegetation invests C to increase N supply in N-limited conditions), and N limitation of decomposition.

We analysed the S3 simulation from the TRENDY protocol, which includes historical changes in atmospheric CO₂, climate, N deposition, N fertilization, and land use from 1851 to 2021 (see Friedlingstein et al., 2022, for a full description of the simulation protocol). Briefly, models were forced with atmospheric CO₂ from Dlugokencky and Tans (2022); the merged monthly Climate Research Unit (CRU), 6-hourly Japanese 55-year Reanalysis (JRA-55) dataset, or the monthly CRU dataset from Harris et al. (2020); N deposition from Hegglin et al. (2016) and Tian et al. (2022); N fertilization from the global N₂O Model Intercomparison Project (NMIP) (Tian et al., 2018); and land use from the LUH2-GCB2022 (Land-Use Harmonization 2) dataset (Chini et al., 2021; Hurtt et al., 2020; Klein Goldewijk et al., 2017a, b). We interpolated outputs from all models to a common resolution of 1^∘ × 1^∘ using bilinear interpolation.

2.2 Terrestrial biosphere model descriptions

The terrestrial biosphere models in the TRENDY ensemble employ a wide variety of assumptions and formulations of N cycling processes, reflecting knowledge gaps and divergent theories (Table 1). Here we describe four fundamental aspects of N cycling for each terrestrial biosphere model: N limitation of vegetation growth, biological N fixation, the response of vegetation to N limitation (i.e., strategies in which vegetation invests C to increase N supply in N-limited conditions), and N limitation of decomposition. These have been identified as important challenges for representing N cycling in terrestrial biosphere models (Meyerholt et al., 2020; Peng et al., 2020; Stocker et al., 2016; Wieder et al., 2015a; Zaehle et al., 2015; Zaehle and Dalmonech, 2011).

Terrestrial biosphere models differ in how N limitation of vegetation growth is represented (Thomas et al., 2015). Some TRENDY models represent flexible C : N stoichiometry and modelled maximum carboxylation rate of photosynthesis (V_cmax) decreases with decreasing leaf N (CABLE-POP, CLASSIC, CLM5.0, LPJ-GUESS, OCNv2, ORCHIDEEv3) following empirical evidence (Walker et al., 2014). Other TRENDY models represent time-invariant C : N stoichiometry and modelled gross primary productivity (GPP) or net primary productivity (NPP) decreases with N limitation (DLEM, ISAM, JSBACH, JULES-ES, and LPX-Bern). Importantly, flexible vs. time-invariant C : N stoichiometry determines terrestrial C storage per unit N.

Biological N fixation is the dominant natural N supply to terrestrial ecosystems (Vitousek et al., 2013). In terrestrial biosphere models, biological N fixation has generally been represented phenomenologically as a function of either NPP or evapotranspiration (ET) (Cleveland et al., 1999). More recently, representations of biological N fixation have been updated such that it is up-regulated in N-limited conditions following empirical evidence (Menge et al., 2015; Vitousek et al., 2013; Zheng et al., 2019). The majority of TRENDY models represent biological N fixation phenomenologically (ISAM, JSBACH, JULES-ES, and LPJ-GUESS). Three TRENDY models (CLASSIC, CLM5.0, and OCNv2) represent biological N fixation mechanistically such that it increases with N limitation of vegetation and has an associated C cost per unit N fixed (Kou-Giesbrecht and Arora, 2022; Lawrence et al., 2019; Meyerholt et al., 2016; Shi et al., 2016; Fisher et al., 2010). These representations separate free-living biological N fixation (via soil microbes, epiphytic microbes, lichens, bryophytes, etc.; Reed et al., 2011) from symbiotic biological N fixation, which is regulated by N limitation of vegetation. DLEM derives biological N fixation as a function of soil temperature, soil moisture, soil C, and soil N. LPX-Bern derives biological N fixation post hoc to simulate a closed N cycle, implicitly including rock N sources (Joos et al., 2020). Finally, CABLE-POP and ORCHIDEEv3 represent biological N fixation as a specified time-invariant input over the historical period. Importantly, representing the regulation of biological N fixation by N limitation not only determines biological N fixation itself but also modulates terrestrial C sequestration: it enables vegetation to increase N uptake in N-limited conditions, reduce N limitation, and sustain terrestrial C sequestration. Some TRENDY models (DLEM, LPJ-GUESS, and OCNv2) also represent increasing C allocation to roots with increasing N limitation (Smith et al., 2014; Zaehle and Friend, 2010) following empirical evidence (Poorter et al., 2012). This enables vegetation to also increase root N uptake in N-limited conditions, reduce N limitation, and sustain terrestrial C sequestration. The response of vegetation to N limitation, which could also include increased C allocation to mycorrhizae (Phillips et al., 2013) (represented in CLM5.0) or increased re-translocation of N during tissue turnover (Du et al., 2020; Han et al., 2013; Kobe et al., 2005) (represented in CLM5.0), is important for determining terrestrial C sequestration.

The decomposition rate is controlled by soil temperature, soil moisture, and N content in litter, where increasing litter C : N ratio decreases the decomposition rate (Cotrufo et al., 2013). Some TRENDY models represent this reduction in decomposition rate with increasing litter C : N ratio (CLM5.0, DLEM, ISAM, JSBACH, JULES-ES, and OCNv2) following empirical evidence.

2.3 Observation-based datasets

We interpolated observation-based datasets to a common resolution of 1^∘ × 1^∘ using bilinear interpolation for comparison against model outputs. To compare model outputs against observation-based datasets we averaged model outputs over 1980–2021, which spans the period in which most measurements were made.

2.3.1 Biological N fixation

A biological N fixation observation-based dataset was derived from Davies-Barnard and Friedlingstein (2020), a global meta-analysis of field measurements of natural biological N fixation (free-living and symbiotic) that scales biome-specific means onto the Collection 5 MODIS Global Land Cover Type International Geosphere-Biosphere Programme (IGBP) product (Friedl et al., 2010). This dataset includes agricultural biological N fixation and assumes that crop biological N fixation rates are equivalent to those of grasses.

The score of LPX-Bern in simulating biological N fixation is not analysed because it implicitly includes rock N sources and is thus not directly comparable to the observation-based dataset.

2.3.2 Vegetation C : N ratio

A vegetation C : N ratio observation-based dataset was derived by scaling biome-specific means for vegetation C : N ratios from the TRY plant trait database (Kattge et al., 2020) onto the Collection 5 MODIS Global Land Cover Type IGBP product (Friedl et al., 2010) and combining it with the remote sensing leaf N content product from Moreno-Martínez et al. (2018). First, we obtained N content per dry mass for leaves, roots, and stems and C content per dry mass for leaves, roots, and stems from the TRY plant trait database. We selected entries that reported species. Second, we obtained the plant functional type (PFT) for each species from the TRY plant trait database. We categorized each PFT into the IGBP land cover types (Table A1 in the Appendix) and then used this to categorize each entry into the IGBP land cover types using species. We averaged across entries in each IGBP land cover type. Third, we divided mean tissue C content per tissue dry mass by mean tissue N content per tissue dry mass for each tissue and for each IGBP land cover type. Fourth, we weighed each tissue by its PFT-specific fraction of total biomass from Poorter et al. (2012) to obtain the total vegetation C : N ratio for each IGBP land cover type. Fifth, we scaled total vegetation C : N ratio and leaf N content per dry mass for each IGBP land cover type to the Collection 5 MODIS Global Land Cover Type IGBP product. Sixth, we multiplied derived total vegetation C : N ratio relative to leaf N content per dry mass by the remote sensing leaf N content per dry mass product (Moreno-Martínez et al., 2018) to obtain a vegetation C : N ratio observation-based dataset.

2.3.3 Soil C : N ratio

A soil C : N ratio observation-based dataset was derived from soil C and soil N products from SoilGrids (Poggio et al., 2021), which provides globally gridded datasets of soil organic C and total soil N at a 250 m × 250 m resolution for six layers up to a depth of 200 cm. These estimates are derived using machine learning methods and soil observations from 240 000 locations across the globe and over 400 environmental covariates. We summed soil C over all layers and soil N over all layers (using the bulk density and depth of each layer) and then obtained the soil C : N ratio.

2.3.4 C cycling variables

In addition to evaluating N cycling variables, we also evaluated the primary C cycling variables: GPP, net biome productivity (NBP), vegetation C (CVEG), soil C (CSOIL), and leaf area index (LAI). These variables have been previously evaluated in detail for the terrestrial biosphere models in the TRENDY ensemble in Seiler et al. (2022). Seiler et al. (2022) give further details on the observation-based datasets used to evaluate the primary C cycling variables. Briefly, we evaluated GPP against MODIS (Zhang et al., 2017), GOSIF (Li and Xiao, 2019), and FLUXCOM (Jung et al., 2020) products. We evaluated NBP against the CAMS (Agustí-Panareda et al., 2019), CarboScope (Rödenbeck et al., 2018), and CT2019 (Jacobson et al., 2020) products. We evaluated CVEG against the GEOCARBON (Avitabile et al., 2016; Santoro et al., 2015), Zhang and Liang (2020), and Huang et al. (2021) products. We evaluated LAI against Advanced Very High Resolution Radiometer (AVHRR; Claverie et al., 2016), Copernicus (Verger et al., 2014), and MODIS (Myneni et al., 2002) products. We evaluated CSOIL against Harmonized World Soil Database (HWSD; Todd-Brown et al., 2013; Wieder, 2014) and SoilGrids (Hengl et al., 2017) products. These observation-based products are globally gridded.

2.4 Model evaluation with the Automated Model Benchmarking R (AMBER) package

The Automated Model Benchmarking R (AMBER) package developed by Seiler et al. (2021) quantifies model performance in reproducing observation-based datasets using a skill score system that is based on ILAMB (Collier et al., 2018). Five scores assess the simulated time mean bias (S_bias), monthly centralized root-mean-square error (S_rmse), seasonality (S_phase), inter-annual variability (S_iav), and spatial distribution (S_dist) in comparison to the observation-based dataset. Scores are dimensionless and range from 0 to 1, where higher values indicate better model performance. The overall score for each variable (S_overall) is

$$\begin{array}{}\text{(1)}& S_{\mathrm{overall}}=\mathrm{mean}\left(S_{\mathrm{bias}},S_{\mathrm{rmse}},S_{\mathrm{phase}},S_{\mathrm{iav}},S_{\mathrm{dist}}\right).\end{array}$$

We calculated the overall score for each C and N cycling variable. Because biological N fixation, vegetation C : N ratio, and soil C : N ratio datasets are representative of the present-day (as a single time point) values, S_rmse, S_phase, and S_iav are not defined and thus do not contribute to S_overall. This also holds for vegetation C and soil C. The calculation of each score is described in detail in Seiler et al. (2022).

2.5 Statistics

We used a Mann–Kendall trend test to assess the existence of a statistically significant trend in the time series over the historical period for simulated C and N cycling variables (Hipel and McLeod, 1994). We conducted two analyses to compare model performance in simulating C cycling vs. N cycling. First, we calculated Spearman's rank correlation coefficient to assess the existence of statistically significant correlations between overall scores, present-day global values, and Kendall's tau. Second, we used a t test or ANOVA (p value < 0.05) to assess the existence of statistically significant differences between overall scores, present-day global values, and Kendall's tau for models with different representations of N limitation of vegetation growth, biological N fixation, vegetation response to N limitation, and N limitation of decomposition (Table 1).

3 Results

3.1 Net biome productivity

Figure 2 shows NBP simulated by the TRENDY ensemble models with coupled C–N cycling (hereafter referred to as the TRENDY-N ensemble). NBP is the difference between the net natural atmosphere–land flux of CO₂ and land use change CO₂ emissions. Positive values of NBP indicate a terrestrial C sink, whereas negative values of NBP indicate a terrestrial C source. All TRENDY-N ensemble models suggest a terrestrial C sink for the present day, agreeing with the global carbon budget constraint from the 2022 Global Carbon Budget with most models within 2 standard deviations of the mean (1.5±0.6 Pg C for 2012–2021) (Fig. 2a). The TRENDY-N ensemble agrees reasonably well with observations globally, agreeing somewhat better with CarboScope and CT2019 than with CAMS (Fig. 2b). However, the latitudinal distributions of the observation-based datasets display weak agreement among themselves with opposing signs in multiple regions due to differences in the inversion models and atmospheric CO₂ measurements used in each dataset (Fig. 2b). The largest differences occur at southern latitudes and at high northern latitudes, and this is in part due to the smaller land area at these latitudes. The regions showing the strongest agreement are at middle to high northern latitudes, where both the TRENDY-N ensemble and observations suggest a terrestrial C sink (Fig. 2b).

Figure 2 Net biome productivity (NBP) simulated by the TRENDY-N ensemble. (a) Global NBP from 1960 to 2021. The boxes indicate the global C budget constraint (difference between fossil fuel CO₂ emissions and the growth rate of atmospheric CO₂ and the uptake of CO₂ by oceans; mean ±2 standard deviation) from the 2022 Global Carbon Budget (Friedlingstein et al., 2022). Thick lines indicate the moving average over 10 years, and thin lines indicate the annual time series. (b) Latitudinal distribution and global mean of NBP (averaged over 1980–2021) in comparison to three datasets (CAMS, Agustí-Panareda et al., 2019; CarboScope, Rödenbeck et al., 2018; and CT2019, Jacobson et al., 2020). The boxplot shows the median, interquartile range (box), and 80 % percentiles (whiskers) of the global mean of NBP.

3.2 Overview of N cycling

Figure 3 shows a schematic of the N cycle alongside the primary N fluxes and C : N ratios of the primary pools simulated by the TRENDY-N ensemble for the present day (averaged over 1980–2021) and observation-based estimates for these variables that have previously been used for model evaluation (Davies-Barnard et al., 2020). Simulated biological N fixation ranged between 20 and 566 Tg N yr⁻¹ (Table 2) in comparison to the observation-based estimate of 88 Tg N yr⁻¹ (52–130 Tg N yr⁻¹). Simulated N₂O emissions ranged between 0.9 and 11.0 Tg N yr⁻¹ (Table 2) in comparison to the observation-based estimate of 10.8 Tg N yr⁻¹ (7.1–16.0 Tg N yr⁻¹) (Tian et al., 2020). Simulated N losses (which include emissions of NH₃, N₂O, NO_x, and N₂, as well as NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ leaching) ranged between 87 and 603 Tg N yr⁻¹ (Table 2) in comparison to the observation-based estimate of 293 Tg N yr⁻¹ (Fowler et al., 2013). The simulated vegetation C : N ratio ranged between 103 and 222 (Table 2) in comparison to the observation-based estimate of 133 (Zechmeister-Boltenstern et al., 2015). The simulated combined litter–soil C : N ratio ranged between 10 and 64 (Table 2) in comparison to the observation-based estimate of 15 (Zechmeister-Boltenstern et al., 2015). Biological N fixation has the largest inter-model spread with a coefficient of variation of 1.06 (Table 2). Figure 4 shows the geographical distribution of the primary N pools and fluxes simulated by the TRENDY-N ensemble for the present day (averaged over 1980–2021), and variation across models is shown in Appendix Fig. A1.

Figure 3 The N cycle and the primary N pools and fluxes simulated by the TRENDY-N ensemble (averaged over 1980–2021). Horizontal black lines indicate observation-based estimates that have previously been used for model evaluation (biological N fixation from Davies-Barnard and Friedlingstein, 2020; vegetation and combined litter-soil C : N ratios from Zechmeister-Boltenstern et al., 2015; N₂O emissions from Tian et al., 2020; and N losses from Fowler et al., 2013). The black box indicates the terrestrial biosphere. N enters the terrestrial biosphere via biological N fixation, N deposition, and N fertilization (entering the organic soil N pool, the inorganic soil N pool (ammonium (NH${}_{\mathrm{4}}^{+}$) or nitrate (NO${}_{\mathrm{3}}^{-})$), or the vegetation N pool). N is transferred from the inorganic soil N pool to the vegetation N pool via N uptake. N is transferred from the vegetation N pool to the litter N pool via N litterfall. N is transferred from the litter N pool to the organic soil N pool via decomposition. N is transferred from the organic soil N pool to the inorganic soil N pool via net N mineralization. N exits the terrestrial biosphere via N loss (which includes N leaching from soils and N₂O, NO_x, NH₃, and N₂ emissions from both soils and land use change). Not all models provide output for each N pool or flux. Note that biological N fixation simulated by LPX-Bern implicitly includes rock N sources.

Figure 4 Geographical distributions of (a) vegetation N, (b) litter N, (c) soil N, (d) biological N fixation, (e) N uptake, (f) net N mineralization, (g) N₂O emissions, and (h) N loss simulated by the TRENDY-N ensemble (averaged across models over 1980–2021). Variation across models is shown in Fig. A1.

Table 2 Global N pools, N fluxes, and C : N ratios simulated by the TRENDY-N ensemble (mean and coefficient of variation across models over 1980–2021).

Figure 5 shows the time series of the change from pre-industrial levels of the primary N pools and fluxes from 1850 to 2021 simulated by the TRENDY-N ensemble. Figure 6 shows the corresponding Kendall's tau, which identifies the existence of a statistically significant trend (Table A2). Over the historical period, some models suggest decreasing vegetation N (6 out of 11 models), whereas other models suggest increasing vegetation N (2 out of 11 models) or no trend in vegetation N (3 out of 11 models). Some models suggest decreasing soil N (7 out of 11 models), whereas other models suggest increasing soil N (4 out of 11 models). Some models suggest increasing biological N fixation (7 out of 11 models), whereas other models suggest decreasing biological N fixation (2 out of 11 models) or no trend in biological N fixation (2 out of 11 models). All models suggest increasing N uptake (10 out of 10 models). Most models suggest increasing net N mineralization rate (9 out of 10 models) or no trend in N mineralization rate (1 out of 10 models). All models suggest increasing N₂O emissions (7 out of 7 models) and increasing N loss (10 out of 10 models).

Figure 5 Time series of the change from the pre-industrial level (averaged over 1850–1870) of (a) vegetation N, (b) litter N, (c) soil N, (d) biological N fixation, (e) N uptake, (f) net N mineralization, (g) N₂O emissions, and (h) N loss simulated by the TRENDY-N ensemble from 1850 to 2021. Figure A5 shows the time series for each N pool and N flux simulated by the TRENDY-N ensemble from 1850 to 2021.

Figure 6 Kendall's tau from the Mann–Kendall test (p value < 0.05) for each N pool and N flux time series simulated by the TRENDY-N ensemble from 1850 to 2021 (Table A2). A positive value (red) indicates an increasing trend and a negative value (blue) indicates a decreasing trend. Gray indicates a statistically insignificant value, and white indicates a missing value.

3.3 Evaluation of biological N fixation, vegetation C : N ratio, and soil C : N ratio

In comparison to the observation-based dataset from Davies-Barnard and Friedlingstein (2020), the TRENDY-N ensemble reproduced global biological N fixation (101.8 Tg N yr⁻¹ vs. 88 Tg N yr⁻¹; Fig. 7a and Table 2) but overestimated low-latitude biological N fixation and underestimated high-latitude biological N fixation in the Northern Hemisphere (Fig. 7b). In comparison to the observation-based dataset from the TRY plant trait database, the TRENDY-N ensemble overestimated the global vegetation C : N ratio (154.5 vs. 102.8; Fig. 7c and Table 2) and overestimated the vegetation C : N ratio across latitudes while capturing its latitudinal pattern (Fig. 7d). In comparison to the observation-based dataset from SoilGrids, the TRENDY-N ensemble overestimated the global soil C : N ratio, simulating a relatively constant soil C : N ratio across latitudes (11.1 vs. 8.8; Fig. 7e and Table 2). The TRENDY-N ensemble was thus unable to capture the latitudinal pattern of the soil C : N ratio (Fig. 7f).

Figure 7 Latitudinal distributions and global means of biological N fixation, vegetation C : N ratio, and soil C : N ratio simulated by the TRENDY-N ensemble (averaged across models over 1980–2021) in comparison to observations. Panels (a, c, e) show the latitudinal distribution of the mean and boxplots show the global mean. Panels (b, d, f) show the latitudinal distribution of the bias. Latitudinal distributions show the mean (black line) and the 50 %, 80 %, and 100 % percentiles across models. Boxplots show the median, interquartile range (box), and 80 % percentiles (whiskers) across models. Observation-based datasets are from Davies-Barnard and Friedlingstein (2020) for biological N fixation, the TRY plant trait database for vegetation C : N ratio, and SoilGrids for soil C : N ratio. LPX-Bern simulations are not shown in (a) or (b). Latitudinal distributions and global means of individual models in the TRENDY-N ensemble are shown in Fig. A6.

The overall score is a metric of model performance in reproducing an observation-based dataset. Overall scores for biological N fixation, vegetation C : N ratio, and soil C : N ratio (0.46, 0.53, and 0.29 averaged across models, respectively) were lower than those for C cycling variables (0.58 averaged across all C cycling variables and across models) (Fig. 8). The mean overall score for vegetation C : N ratio across models (0.53) was lower than the mean overall scores for vegetation C across models (which ranged from 0.61 to 0.69 depending on the observation-based dataset used to derive the score). Similarly, the mean overall score for soil C : N ratio across models (0.29) was lower than the mean overall scores for soil C across models (which ranged from 0.39 to 0.53 depending on the observation-based dataset used to derive the score).

Figure 8 Overall scores of the TRENDY-N ensemble in simulating C and N cycling variables: gross primary productivity (GPP), net biome productivity (NBP), vegetation C (CVEG), soil C (CSOIL), leaf area index (LAI), biological N fixation (FBNF), vegetation C : N ratio (CNVEG), and soil C : N ratio (CNSOIL). Abbreviations of the observation-based datasets are described in Sect. 2 and in Seiler et al. (2022).

For N cycling variables, the overall score is composed of the time mean bias score (which assesses the difference between the time mean of model simulations and the time mean of the observation-based dataset) and the spatial distribution score (which assesses the ability of the model to reproduce the spatial pattern of the observation-based dataset) (Collier et al., 2018; Seiler et al., 2022). For biological N fixation, the time mean bias score averaged across models was 0.50 and the mean spatial distribution score across models was 0.41 (Table A3). For the vegetation C : N ratio, the time mean bias averaged score across models was 0.46 and the mean spatial distribution score across models was 0.59 (Table A3). For the soil C : N ratio, the time mean bias score averaged across models was 0.39 and the mean spatial distribution score across models was 0.19 (Table A3).

Note that for C fluxes the overall score is composed of not only the time mean bias score and the spatial distribution score but also the monthly centralized root-mean-square-error score (which assesses the ability of the model to reproduce the time series of the observation-based dataset), the seasonality score (which assesses the ability of the model to reproduce the seasonality of the observation-based dataset), and the inter-annual variability score (which assesses the ability of the model to reproduce the inter-annual variability of the observation-based dataset) because observation-based datasets of C fluxes are available over time, whereas observation-based datasets of C pools and all N cycling variables are representative of the present day (as a single time point).

3.4 Model performance for C cycling vs. N cycling

There were no statistically significant correlations between the overall score of NBP (as well as other primary C variables) and the overall scores of the primary N variables across the TRENDY-N ensemble (Fig. A2). Furthermore, there were no statistically significant correlations between the present-day global value of NBP and the present-day global values of the primary N variables across the TRENDY-N ensemble (Fig. A3). Finally, there were no statistically significant correlations between Kendall's tau of NBP and Kendall's tau of the primary N variables across the TRENDY-N ensemble (Fig. A4).

3.5 Model performance for different representations of N cycling processes

There were no statistically significant differences in overall scores between models with different representations of N limitation of vegetation growth (decreasing V_cmax and flexible C : N stoichiometry vs. decreasing NPP), different representations of biological N fixation (function of N limitation of vegetation growth vs. function of NPP or ET vs. time invariant), different representations of the response of vegetation to N limitation (dynamic vs. static), or different representations of N limitation of decomposition (function of soil N vs. N invariant) (Table A4). However, models that represented decomposition as a function of soil N had a significantly higher NBP score (for CT2019) than models that represented decomposition as N invariant. Similarly, there were no statistically significant differences between present-day global values or Kendall's tau of primary C and N pools and fluxes between models with different representations of N limitation of vegetation growth, biological N fixation, vegetation response to N limitation, and N limitation of decomposition (Tables A5 and A6). This is likely in part due to the low number of models and the confounding influence of other process representations.

4 Discussion

4.1 Evaluation of N cycling in terrestrial biosphere models

Despite the ability of all TRENDY-N models to simulate the historical terrestrial C sink in line with observations (Fig. 2), there is substantial variation in simulated N cycling processes by the models. The magnitudes of N pools and fluxes differ considerably between models (Figs. 3 and A1). Additionally, the historical trajectories of these N pools and fluxes differ between models: some models simulate increasing vegetation N and soil N, whereas others simulate decreasing vegetation N and soil N between 1850 and 2021 (Figs. 5 and 6). These trajectories are the result of a host of interacting global change drivers (CO₂ fertilization, intensifying N deposition, rising temperature and varying precipitation, and land use change and associated N fertilization regimes) whose effects are challenging to disentangle without additional simulations. For example, while intensifying N deposition and N fertilizer use could drive increasing soil N and N uptake, land use change could increase N losses from both vegetation N and soil N. Most models suggest increasing biological N fixation between 1850 and 2021. This occurs as a result of either increasing vegetation biomass or the up-regulation of biological N fixation due to N limitation imposed by CO₂ fertilization or a combination thereof, depending on the representation of biological N fixation in a given model (Table 1). This follows observations that suggest that biological N fixation is stimulated by CO₂ fertilization (Zheng et al., 2020; Liang et al., 2016), although its mechanism (i.e., up-regulated biological N fixation in N-limited conditions) may not be captured. Similarly, most models also suggest increasing N uptake between 1850 and 2021. This also occurs as a result of increasing vegetation biomass; increasing soil N from intensifying N deposition and N fertilizer use; or increasing biological N fixation, mycorrhizae, and root allocation due to N limitation imposed by CO₂ fertilization, which is again dependent on the representation of the vegetation response to N limitation in a given model (Table 1). Most models suggest increasing net N mineralization rate between 1850 and 2021, likely due to rising temperature following observations (Liu et al., 2017). Most models suggest increasing N₂O emissions (and N losses) between 1850 and 2021, likely due to rising temperatures and intensifying N deposition and N fertilizer use following observations (Tian et al., 2020).

We focused on three key N cycling processes for evaluation: biological N fixation, vegetation C : N ratio, and soil C : N ratio. These three key N cycling processes have important implications for projecting the future terrestrial C sink. Biological N fixation is the dominant natural N supply to terrestrial ecosystems and allows vegetation to increase N uptake in N-limited conditions, reduce N limitation, and thus sustain terrestrial C sequestration, such as in response to N limitation imposed by CO₂ fertilization (Zheng et al., 2020; Liang et al., 2016). Vegetation and soil C : N ratios reflect assimilated C per unit N and thus terrestrial C sequestration. They can potentially vary, such as in response to high photosynthesis rates relative to N uptake rates driven by CO₂ fertilization (Elser et al., 2010). Overall scores of N cycling variables, which quantify model performance in reproducing an observation-based dataset, are lower than overall scores of corresponding C cycling variables, suggesting that models could be less capable of capturing N cycling processes than C cycling processes (Fig. 8). However, this could also be due to the significant uncertainty associated with measurements of N cycling processes as discussed below.

The TRENDY-N ensemble reproduced global observation-based biological N fixation but tended to overestimate low-latitude biological N fixation and underestimate high-latitude biological N fixation (Fig. 7a, b). This is likely because most models represented biological N fixation phenomenologically as a function of a measure of vegetation activity (either NPP or ET). Since there is higher vegetation activity at low latitudes than at high latitudes, these models thus represent higher biological N fixation at low latitudes than at high latitudes. However, because biological N fixation is down-regulated in non-N-limited conditions, it is often down-regulated at low latitudes, which are generally not (or at least less) N-limited in nature (Barron et al., 2011; Batterman et al., 2013; Sullivan et al., 2014). While CLASSIC, CLM5.0, and OCNv2 can represent the down-regulation of biological N fixation in non-N-limited conditions, they still simulate high low-latitude biological N fixation. This suggests that the strength of regulation of biological N fixation could be insufficient and/or that there could be unaccounted N sources at low latitudes. For example, rock N weathering could be a significant N source to terrestrial ecosystems. Some estimates have suggested that rock N weathering could be as high as 11–18 Tg N yr⁻¹ globally (Houlton et al., 2018) but is not explicitly represented in the TRENDY-N ensemble (with the exception of LPX-Bern, which calculates all external N sources post hoc to simulate a closed N cycle, thereby implicitly including rock N sources). The discrepancy between modelled and observed biological N fixation could also be due to uncertainty in the observation-based dataset given the difficulties associated with measuring biological N fixation (Soper et al., 2021). Ecological theory (Hedin et al., 2009) has suggested that natural biological N fixation should be higher at low latitudes given large N losses, in contrast to the observation-based dataset from Davies-Barnard and Friedlingstein (2020). Furthermore, the observation-based dataset from Davies-Barnard and Friedlingstein (2020) did not explicitly account for agricultural biological N fixation but rather assumed that crop biological N fixation rates are equivalent to those of grasses, although they are likely to be much greater (Peoples et al., 2021; Herridge et al., 2022).

The TRENDY-N ensemble overestimated the global observation-based vegetation C : N ratio but reproduced its latitudinal pattern (as also indicated by its higher spatial distribution score) (Fig. 7c, d). This is because most models represent different plant functional types (e.g., evergreen needleleaf trees, deciduous broadleaf trees, evergreen broadleaf trees) with different tissue C : N ratios (which can be either flexible within a constrained range or time invariant). These plant functional types are geographically distributed according to similar land cover products. The TRENDY-N ensemble overestimated global observation-based soil C : N ratio and failed to reproduce its latitudinal pattern (as also indicated by its lower spatial distribution score) (Fig. 7e, f). In particular, models failed to reproduce the peak at the Equator and the peak at approximately −30^∘ S, corresponding to tropical forests and deserts, respectively. This is because most models represent a constant soil C : N ratio (both temporally and spatially) and are thus unable to capture the spatial variability in the soil C : N ratio. Improving the representation of soil N is an important future direction for terrestrial biosphere model development given the essential feedbacks between soil N and soil C.

4.2 Disconnect between C and N cycling in terrestrial biosphere models

The importance of N limitation of terrestrial C sequestration is empirically established (Elser et al., 2007; LeBauer and Treseder, 2008; Wright et al., 2018). It has already influenced the historical terrestrial C sink (S. Wang et al., 2020), and it is expected to be especially important under future CO₂ fertilization and global change (Terrer et al., 2019). While all TRENDY-N models simulate the historical terrestrial C sink in line with observations (and are no different from TRENDY models without a representation N cycling; Seiler et al., 2022), our results suggest a disconnect between C and N cycling in these models. First, the models exhibit a wide spread across simulated N pools and fluxes. Second, there are no significant correlations between model performance in simulating N cycling and model performance in simulating C cycling. Third, there are no statistically significant differences between models with different representations of fundamental N cycling processes (N limitation of vegetation growth, biological N fixation, the response of vegetation to N limitation, and N limitation of decomposition).

Overall, our results suggest that the underlying N cycling processes that regulate terrestrial C sequestration operate differently across models and may not be fully captured given that models are calibrated to C cycling. The spread across models suggests that approaches to represent N cycling processes vary among models and that there is no clear consensus yet on what the best approaches are. Studies have explored the validity of different representations of N cycling processes within a single model, suggesting that alternative representations of a biological N fixation, ecosystem C : N stoichiometry, and ecosystem N losses lead to substantial differences in simulated C cycling (Kou-Giesbrecht and Arora, 2022; Meyerholt et al., 2020; Peng et al., 2020; Wieder et al., 2015a). This disconnect between C and N cycling will become particularly consequential for projecting the terrestrial C sink under future global change, which is likely to modify the C–N balance through N limitation of CO₂ fertilization and intensifying N deposition (among other effects of global change).

4.3 Future directions

Evaluating N cycling in terrestrial biosphere models is severely restricted by the lack of available observations of N cycling. N cycling processes are notoriously difficult to measure, such as biological N fixation (Soper et al., 2021) and gaseous N losses (Barton et al., 2015). In the past, N cycling has commonly been evaluated by comparison to estimates of global N pools and fluxes derived from a small number of observations that have been scaled up or averaged to yield a value with wide confidence intervals (Davies-Barnard et al., 2020). Not only are these global totals highly uncertain, but they also do not allow for the analysis of spatial patterns. Here, we present an improved framework to evaluate three key N cycling processes – biological N fixation, vegetation C : N ratio, and soil C : N ratio – in terrestrial biosphere models. However, these globally gridded observation-based datasets are also uncertain, given uncertainty in the estimates of tissue C : N ratios for different plant functional types, tissue fraction of total biomass (especially those of roots and wood, which had a lower number of measurements in comparison to that of leaves), and the measurements and models used to derive soil N (Batjes et al., 2020). More observations of these N cycling processes are necessary to reduce uncertainty. Temporally explicit measurements are important for assessing intra-annual and inter-annual variability. Leveraging advances in remote sensing (Knyazikhin et al., 2013; Townsend et al., 2013; Cawse-Nicholson et al., 2021) and incorporating N cycling process measurements into research networks such as FLUXNET (Vicca et al., 2018) are essential. Multiple observation-based datasets taken from different sources and derived via different methodologies of a given N cycling process are necessary to evaluate observational uncertainty (Seiler et al., 2021). Global observations of other important N cycling processes (such as N mineralization and N losses) are necessary to fully evaluate N cycling in terrestrial biosphere models. Additionally, hindcast simulations of the transition from the Last Glacial Maximum to the pre-industrial period can be used in combination with proxy-based reconstructions of past N₂O emissions (Fischer et al., 2019) and C stocks (Jeltsch-Thömmes et al., 2019) for model evaluation and can serve as a constraint for terrestrial biosphere models (Joos et al., 2020).

Modelled experimental manipulations (such as CO₂ fertilization or N fertilization experiments) are imperative to evaluate model formulations of the underlying mechanisms of C–N cycling interactions (Medlyn et al., 2015; Wieder et al., 2019; Zaehle et al., 2014). Derived nutrient limitation products (Fisher et al., 2012) can also be applied to evaluate present-day nutrient cycling when phosphorus (P) is accounted for (Braghiere et al., 2022). Evaluating the ability of models to simulate present-day N cycling processes, as we did here, is only one method of assessing their ability to simulate N limitation of terrestrial C sequestration. A robust test of the simulated response to CO₂ fertilization and N fertilization across models would be ideal for evaluating the ability of models to represent the regulation of C cycling by N cycling under global change and thus their ability to realistically simulate the future terrestrial C sink.

While some of the models in the TRENDY-N ensemble have the capability of representing coupled C, N, and P cycling (Goll et al., 2012; Nakhavali et al., 2022; Sun et al., 2021; Wang et al., 2010; Z. Wang et al., 2020; Yang et al., 2014), P cycling was not active in the model simulations in the GCP 2022. P limitation could be important for limiting terrestrial C sequestration, especially in low-latitude forests (Elser et al., 2007; Terrer et al., 2019; Wieder et al., 2015b). As more models incorporate coupled C–N–P cycling (Reed et al., 2015; Braghiere et al., 2022), observation-based datasets of P will also be necessary for model evaluation.

5 Conclusions

Because the TRENDY-N ensemble overestimated both vegetation and soil C : N ratios, it is possible that models could overestimate assimilated C per unit N and thus future terrestrial C sequestration under CO₂ fertilization. Alongside discrepancies in biological N fixation, this could lead to biases in projections of the future terrestrial C sink by the TRENDY-N ensemble. Not to mention there are several other terrestrial biosphere models in the TRENDY ensemble that do not represent coupled C–N cycling. While the models are capable of reproducing the current terrestrial C sink, the spread across the models in simulating N cycling suggests that C–N interactions operate differently across models and may not be fully captured given that models are calibrated to C cycling. However, these C–N interactions are critical for projecting the terrestrial C sink under global change in the future.

Appendix A

Table A1 IGBP land cover type, corresponding TRY plant trait database PFT, tissue C : N ratios from the TRY plant trait database (Kattge et al., 2020), tissue fractions (Poorter et al., 2012), and calculated total C : N ratio.

^a Value from evergreen needleleaf forest. ^b Average of evergreen needleleaf forest, evergreen broadleaf forest, and deciduous broadleaf forest. ^c Value from grasslands. ^d Value from croplands.

Table A2 Kendall's tau from the Mann–Kendall test (p value < 0.05) for each N pool and N flux time series simulated by the TRENDY-N ensemble from 1850 to 2021. NS indicates that Kendall's tau is not significant. NA indicates that the variable was not reported by the model.

Table A3 Time mean bias score (S_bias), spatial distribution score (S_dist), and overall score (S_overall) of the TRENDY-N ensemble in simulating biological N fixation, vegetation C : N ratio, and soil C : N ratio. NA indicates that this variable was not evaluated.

Table A4 Overall scores of biological N fixation, vegetation C : N ratio, soil C : N ratio, and NBP averaged across TRENDY-N ensemble models with different representations of key N cycling processes (N limitation of vegetation growth, biological N fixation, vegetation response to N limitation, and N limitation of decomposition; see Table 1). The p values are from t tests and ANOVAs assessing differences between these representations of key N cycling processes.

Table A5 Present-day global values of biological N fixation, vegetation C : N ratio, and soil C : N ratio averaged across TRENDY-N ensemble models with different representations of key N cycling processes (N limitation of vegetation growth, biological N fixation, vegetation response to N limitation, and N limitation of decomposition; see Table 1). The p values are from t tests and ANOVAs assessing differences between these representations of key N cycling processes.

Table A6 Kendall's tau from the Mann–Kendall test (p value < 0.05) for biological N fixation, vegetation C : N ratio, and soil C : N ratio averaged across TRENDY-N ensemble models with different representations of key N cycling processes (N limitation of vegetation growth, biological N fixation, vegetation response to N limitation, and N limitation of decomposition; see Table 1). The p values are from t tests and ANOVAs assessing differences between these representations of key N cycling processes.

Figure A1 Geographical distributions of variation in (a) vegetation N, (b) litter N, (c) soil N, (d) biological N fixation, (e) N uptake, (f) net N mineralization, (g) N₂O emissions, and (h) N loss simulated by the TRENDY-N ensemble (across models over 1980–2021).

Figure A2 Correlations between overall scores of primary C and N pools and fluxes across TRENDY-N ensemble models: gross primary productivity (GPP), net biome productivity (NBP), vegetation C (CVEG), soil C (CSOIL), leaf area index (LAI), biological N fixation (FBNF), vegetation C : N ratio (CNVEG), and soil C : N ratio (CNSOIL). Abbreviations of the observation-based datasets are described in Sect. 2 and in Seiler et al. (2022). Spearman's rank correlation coefficient is shown for statistically significant correlations (p value < 0.05).

Figure A3 Correlations between present-day global values (averaged over 1980–2021) of primary C and N pools and fluxes across TRENDY-N ensemble models: vegetation C (CVEG), litter C (CLITTER), soil C (CSOIL), net biome productivity (NBP), gross primary productivity (GPP), autotrophic respiration (RA), heterotrophic respiration (RH), leaf area index (LAI), vegetation N (NVEG), litter N (NLITTER), soil N (NSOIL), biological N fixation (FBNF), N uptake (NUP), net N mineralization (NETNMIN), N₂O emissions (N₂O), N loss (NLOSS), vegetation C : N ratio (CNVEG), and soil C : N ratio (CNSOIL). Spearman's rank correlation coefficient is shown for statistically significant correlations (p value < 0.05).

Figure A4 Correlations between Kendall's tau of primary C and N pools and fluxes across TRENDY-N ensemble models: vegetation C (CVEG), litter C (CLITTER), soil C (CSOIL), net biome productivity (NBP), gross primary productivity (GPP), autotrophic respiration (RA), heterotrophic respiration (RH), leaf area index (LAI), vegetation N (NVEG), litter N (NLITTER), soil N (NSOIL), biological N fixation (FBNF), N uptake (NUP), net N mineralization (NETNMIN), N₂O emissions (N₂O), N loss (NLOSS), vegetation C : N ratio (CNVEG), and soil C : N ratio (CNSOIL). Spearman's rank correlation coefficient is shown for statistically significant correlations (p value < 0.05).

Figure A5 Time series of (a) vegetation N, (b) litter N, (c) soil N, (d) biological N fixation, (e) N uptake, (f) net N mineralization, (g) N₂O emissions, and (h) N loss simulated by the TRENDY-N ensemble from 1850 to 2021.

Figure A6 Latitudinal distributions and global means of (a) biological N fixation, (b) vegetation C : N ratio, and (c) soil C : N ratio simulated by the TRENDY-N ensemble (averaged across models over 1980–2021) in comparison to observation-based datasets from Davies-Barnard and Friedlingstein (2020) for biological N fixation, the TRY plant trait database for vegetation C : N ratio, and SoilGrids for soil C : N ratio. Boxplots show the median, interquartile range (box), and 80 % percentiles (whiskers) of the global mean.