Biological nitrogen fixation of natural and agricultural vegetation simulated with LPJmL 5.7.9

Biological nitrogen fixation (BNF) by symbiotic and free-living bacteria is an important source of plant-available nitrogen (N) in terrestrial ecosystems supporting carbon (C) sequestration and food production worldwide. Dynamic global vegetation models (DGVMs) are frequently used to assess the N and C cycles under dynamic land use and climate. BNF plays an important role in the components of both these cycles, making a robust representation of the processes and variables that BNF depends on important to reduce uncertainty within the C and N cycles and improve the ability of DGVMs to project future ecosystem productivity, vegetation patterns or the land C sink. Still, BNF is often modelled as a function of net primary productivity or evapotranspiration, and the actual drivers are neglected. We implemented plant-functional-type-specific limitations for BNF dependent on soil temperature and soil water content, as well as a cost of BNF, in the Lund–Potsdam–Jena managed Land (LPJmL) DGVM and compared the new (“C-costly”) against the previous (“Original”) approach and data from the scientific literature. For our comparison, we simulated a potential natural vegetation scenario and one including anthropogenic land use for the period from 1901 to 2016 for which we evaluate BNF and legume crop yields. Our results show stronger agreement with BNF observations for the C-costly than the Original approach for natural vegetation and agricultural areas. The C-costly approach reduced the overestimation of BNF, especially in hot spots of legume crop production. Despite the reduced BNF in the C-costly approach, yields of legume crops were similar to the Original approach. While the net C and N balances were similar between the two approaches, the reduced BNF in the C-costly approach results in a slight underestimation of N losses from leaching, emissions and harvest compared to the values in the literature, supporting further investigation of the underlying reasons, such as processes represented in DGVMs and scenario assumptions. While we see the potential for further model development, for example, to separate symbiotic and free-living BNF, the C-costly approach is a major improvement over the simple Original approach because of the separate representation of important drivers and limiting factors of BNF, and the C-costly approach also improves the ability of LPJmL to project future C and N cycle dynamics.

1 Introduction

Biological nitrogen fixation (BNF) is an important source of plant-available nitrogen (N) in terrestrial ecosystems (Galloway et al., 1995). It can be separated into symbiotic (Granhall, 1981) and free-living (Reed et al., 2011) BNF, which account for the total BNF with different shares in different ecosystems (Davies-Barnard and Friedlingstein, 2020b). In natural terrestrial ecosystems, N deposition, N fixation through lightning, and BNF are the only processes that introduce additional reactive N into the system (Yu and Zhuang, 2020). In agricultural systems, increased N inputs are – together with extensive manure recycling – a major source of nitrous oxide (N₂O) and ammonium (${{\mathrm{NH}}_{\mathrm{4}}}^{+}$) emissions (Reay et al., 2012; Tian et al., 2020) and nitrate (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$) pollution (Moss, 2007). These inputs result from increased BNF and the deposition of additional anthropogenic N inputs, which originate mainly from synthetic fertiliser application (Lu and Tian, 2017). Promoting N-fixing crops such as forage and grain legumes for usage as green manure has been discussed (Becker et al., 1995; Fageria, 2007; Northup and Rao, 2016) to reduce N losses from nitrification, volatilisation, denitrification and leaching on agricultural land. Generally, symbiotic BNF, as well as free-living BNF, can be important for plant growth in N-limited ecosystems, and this supports carbon (C) sequestration and food production across the globe.

Briefly, BNF describes the transformation of atmospheric N₂ to ammonia (${{\mathrm{NH}}_{\mathrm{4}}}^{+}$) by a variety of soil microorganisms providing a source of mineral N for plants at the expense of C (Yu and Zhuang, 2020). The underlying mechanisms of BNF, as well as its role within the C and N cycles and for ecosystem productivity, have been described in detail in multiple studies (e.g. Yu and Zhuang, 2020; Davies-Barnard and Friedlingstein, 2020b; Cleveland et al., 1999). Here, we focus on the representation of BNF in the Lund–Potsdam–Jena managed Land (LPJmL) dynamic global vegetation model (DGVM) (Schaphoff et al., 2018b; von Bloh et al., 2018; Lutz et al., 2019; Herzfeld et al., 2021; Porwollik et al., 2022; Heinke et al., 2023). We do not distinguish between symbiotic and free-living BNF throughout this study but only consider the total BNF as the sum of both forms.

DGVMs such as LPJmL can be used to assess the role of BNF for the productivity of natural and agricultural ecosystems and its effects on the N and C cycles under dynamic land use and climate. A solid representation of the processes behind BNF is important to reduce uncertainty and improve the model results of DGVMs, which are frequently used in impact assessments and to inform policy-makers. A variety of approaches of different complexity to model BNF have been developed. A key difference between approaches is the selection of variables that control BNF and the accounting of the C cost of BNF. For example, Cleveland et al. (1999) use actual evapotranspiration as a single explanatory variable, while Yu and Zhuang (2020) consider soil temperature, soil water content, soil mineral N and soil C content. Both of these approaches do not consider the cost of BNF neglecting the reduced C assimilation (Cleveland et al., 1999; Yu and Zhuang, 2020), while others explicitly consider a cost per amount of N fixed and a maximum amount of C that can be invested in BNF (e.g. Ma et al., 2022). Even more complex approaches consider the different pathways of N uptake that are associated with a cost (active N uptake, retranslocation and BNF) and optimise for the minimum cost (e.g. Fisher et al., 2010). Depending on the considered variables, the simulated BNF and how it is affected by climate change may strongly differ, which in turn can have strong effects on the simulated C and N fluxes and pools.

A comparison to data published by Davies-Barnard and Friedlingstein (2020a) suggests that the approach that was implemented in LPJmL (von Bloh et al., 2018) based on Cleveland et al. (1999) – in the following defined as the “Original” approach – overestimates global BNF. In addition, we identified several shortcomings of the Original approach in LPJmL. In the Original approach, BNF is a function of actual evapotranspiration, which leads to an overestimation of BNF in moist but not necessarily N-limited ecosystems and an underestimation in dry but N-limited ecosystems. In this simplified implementation, BNF is not constrained by the availability of reactive forms of N, and additional N is fixed even if the reactive soil N is sufficient to fulfil the N demand, which potentially leads to an overestimation of the ammonia pool and N losses. For cultivated grain legumes, the approach assumes no limitation of BNF at all but simply supplies all N requested by the plant that cannot be fulfilled through N uptake from mineral N pools in the soil. This leads to an overestimation of cropland BNF. In order to overcome these deficiencies, here we describe a revision of the Original approach in LPJmL with a more complex approach, referred to as “C-costly” approach in the following. The C-costly approach is inspired by Ma et al. (2022) and Yu and Zhuang (2020) and introduces plant-functional-type (PFT)-specific limitations for BNF dependent on soil temperature and soil water content, as well as a C cost of BNF. In the following, we present the C-costly BNF approach and evaluate its performance against global and site-specific data. We discuss the differences between the Original and the C-costly BNF approach for the N cycle and plant productivity.

2 Methods

2.1 Model description

LPJmL is a dynamic global vegetation model (DGVM) with the full terrestrial hydrology and explicit representation of agricultural management systems for cropland and pastures. We have implemented the BNF module in the most recent development branch, which is based on a consolidated version of the carbon-only model (LPJmL4, Schaphoff et al., 2018b, a), the N cycle (LPJmL5, von Bloh et al., 2018), tillage (Lutz et al., 2019), manure (Herzfeld et al., 2021), cover crop (Porwollik et al., 2022) and grazing management (Heinke et al., 2023) modules. There have been further model improvements that have not been described in publications elsewhere, including improved online coupling options with other models such as IMAGE (Müller et al., 2016) or copan:CORE (Donges et al., 2020). For a better representation of crops that are not explicitly represented (referred to as “others”), these are no longer assumed to be identical to managed grassland (Bondeau et al., 2007) but can be simulated as separate stands with distinct management inputs (e.g. fertiliser amounts).

The original spin-up protocol for LPJmL4, described in Schaphoff et al. (2013), was modified to account for the interaction between soils and plants through N supply in LPJmL5. The principal technique to accelerate the spin-up by calculating the equilibrium soil C stocks from litter decomposition (i.e. the flux of C into the soil C pools) and soil C turnover rates (or residence time) remains the same as in Schaphoff et al. (2013). However, the original code was refactored to improve the accuracy of estimates of equilibrium stocks and to apply the technique to soil C and N pools simultaneously.

In LPJmL5, an adjustment of N pools can lead to a change in plant productivity through a change in N supply from mineralisation. To account for this feedback, the C and N stock adjustments need to be repeated multiple times until the soil and the vegetation reach equilibria. The revised spin-up procedure starts with an initial period of 300 years, during which vegetation is allowed to establish. This is followed by a 2400-year period, during which soil C and N pools are updated every 15 years based on the litter decomposition and soil pool turnover rates of the preceding 10 years. This long period with repeated adjustment (160 times) of C and N pools is required to reach an equilibrium in regions with very low turnover rates (e.g. in the boreal zone). To reduce the effect of inter-annual variability on estimates of equilibrium stocks, a final adjustment is applied after 300 simulation years, using the litter decomposition and soil pool turnover rates over that period. Finally, the model is allowed to adjust to the new C and N stocks for another 500 simulation years.

To assess the effectiveness of the spin-up procedure, we conducted a 1000-year model run under the same conditions as during the spin-up period (i.e. stable pre-industrial atmospheric CO₂ concentration, atmospheric N deposition and climate) for which we present results in Appendix C.

Further changes to the code since the last published version (see Porwollik et al., 2022) include various bug fixes concerning fertiliser and manure application, data output, environmental flow requirements (Jägermeyr et al., 2017), soil temperature (Schaphoff et al., 2013) and bioenergy plantations (Beringer et al., 2011). The latest code changes are now also documented in a CHANGELOG.md file as part of the code repository (Wirth et al., 2024).

2.2 BNF-relevant nitrogen cycle components in LPJmL

While we refer to von Bloh et al. (2018) for a detailed description and evaluation of the N cycle in LPJmL, we briefly describe the main processes that determine N deficit – which is the prerequisite for N fixation in the C-costly approach – and the Original approach here and provide the full equations in Appendix A. An N deficit is defined as the difference between the plant N demand (Eq. A4) and the active and passive N-uptake (Eq. A5) and labile-N reserves (Eq. A11).

$$\begin{array}{}\text{(1)}& N_{\text{deficit},t}=N_{\text{demand},t}-(N_{\text{uptake},t}+N_{\text{labile},t})\end{array}$$

The N demand accounts for N required to produce RuBisCo, depending on the maximum carboxylation capacity and the leaf area index (LAI) of the respective PFT (Eq. A1, first summand) and the structural N demand, depending on the current N content of the different plant compartments (Eq. A1, second summand, and Eq. A4). N reserves are included using a PFT-specific parameter (Eq. A4).

The N uptake is calculated as a combination of passive and active N uptake from the soil and is a function of the potential N uptake of the root system (Eq. A5), which is reduced to account for soil mineral N availability (Eq. A8), soil temperature (Eq. A9) and plant N starvation (Eq. A10). Labile N reserves represent the N currently available from past N uptake, BNF or retranslocation (Eq. A11).

In the Original approach, BNF was calculated from the 20-year average of annual evapotranspiration (etp) for tree and herbaceous PFTs, following the function from Cleveland et al. (1999):

$$\begin{array}{}\text{(2)}& \text{BNF}=\left\{\begin{array}{ll}max(\mathrm{0},(\mathrm{0.0234}\cdot \text{etp}-\mathrm{0.172})/\mathrm{10}/\mathrm{365})& \\ & \mathrm{\text{if }}{C}_{\text{root}}>\mathrm{20}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}\\ \mathrm{0}& \mathrm{\text{otherwise}}.\end{array}\right.\end{array}$$

The resulting BNF is added to the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool of the first soil layer. For crop PFTs, BNF equals N_deficit and is directly added to N_labile.

2.3 The C-costly approach

A key feature is the connection of BNF to an associated cost represented as a reduction in the net primary production (NPP). The C-costly approach calculates actual BNF (N_fix) from the potential BNF (N_fix,pot) using several reduction factors. First, the N fixation rate for the environmental conditions N_fix,env is calculated from N_fix,pot for the first two soil layers, $l=\mathrm{1},\mathrm{2}$, accounting for reductions by dimensionless soil temperature and soil water content (SWC) limitation functions (f_T, f_W) and the root distribution rootdist in the interval [0,1] (Ma et al., 2022):

$$\begin{array}{}\text{(3)}& N_{\text{fix,env}}=\sum _{l=\mathrm{1}}^{\mathrm{2}}{N}_{\text{fix,pot}}\cdot f_T\left(T_{\text{soil},l}\right)\cdot f_{\mathrm{W}}\left({\text{SWC}}_l\right)\cdot {\text{rootdist}}_l.\end{array}$$

The soil temperature limitation is increasing linearly outside the optimal temperature interval, [T_opt,low,T_opt,high] (Eq. 4; Fig. 1a), and it prohibits BNF if outside the tolerable temperature interval, [T_min,T_max], while the soil water limitation is linearly dependent on the relative soil water content, SWC (Eq. 5; Fig. 1b).

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{f}_T& {\displaystyle }\left(T_{\text{soil}}\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{0},& \text{if }{T}_{\text{soil}}<T_{\text{min}}\text{ or }{T}_{\text{soil}}>T_{\text{max}}\\ \frac{T_{\text{soil}}-T_{\text{min}}}{T_{\text{opt,low}}-T_{\text{min}}},& \text{if }{T}_{\text{min}}\le T_{\text{soil}}<T_{\text{opt,low}}\\ \mathrm{1},& \text{if }{T}_{\text{opt,low}}\le T_{\text{soil}}\le T_{\text{opt,high}}\\ \frac{T_{\text{max}}-T_{\text{soil}}}{T_{\text{max}}-T_{\text{opt,high}}},& \text{if }{T}_{\text{opt,high}}<T_{\text{soil}}\le T_{\text{max}}.\end{array}\right.\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{f}_{\mathrm{W}}& {\displaystyle }\left(\text{SWC}\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{0},& \text{if }\text{SWC}\le {\text{SWC}}_{\text{low}}\\ {\mathit{\phi }}_{\mathrm{1}}+\text{SWC}\cdot {\mathit{\phi }}_{\mathrm{2}},& \text{if }{\text{SWC}}_{\text{low}}<\text{SWC}<{\text{SWC}}_{\text{high}}\\ \mathrm{1},& \text{if }\text{SWC}\ge {\text{SWC}}_{\text{high}}.\end{array}\right.\end{array}\end{array}$$

The root distribution is calculated as in Eq. (A7). Since only the fraction of roots in the first two soil layers is used for BNF, shallow root profiles lead to a higher BNF compared to deep root profiles. N_fix,pot, T_min, T_opt,low, T_opt,high, T_max, SWC_low, SWC_high, φ₁ and φ₂ are PFT-specific parameters (Table 1), and their values are adopted from Yu and Zhuang (2020) for the natural vegetation PFTs and from Ma et al. (2022) for soybean and pulses.

Figure 1 Dimensionless temperature limitation function f_T(T) (a) and soil water limitation function f_W(SWC) (b).

Table 1 BNF-related PFT-specific parameter values for the tropical broadleaved evergreen tree (TrBE), tropical broadleaved raingreen tree (TrBR), temperate needleleaved evergreen tree (TeNE), temperate broadleaved evergreen tree (TeBE), temperate broadleaved summergreen tree (TeBS), boreal needleleaved evergreen tree (BoNE), boreal broadleaved summergreen tree (BoBS), boreal needleleaved summergreen tree (BoNS), tropical herbaceous (TrH), temperate herbaceous (TeH), polar herbaceous (PoH), soybean and pulses.

If N_fix,env exceeds the amount of N missing to fulfil the N demand of the current day (the N deficit N_deficit), the N fixation is reduced as follows:

$$\begin{array}{}\text{(6)}& N_{\text{fix,need}}=min(N_{\text{deficit}},N_{\text{fix,env}}).\end{array}$$

Finally, if the cost of the N fixation exceeds the NPP available for BNF, then the N fixation is further reduced to match the maximum amount that can be fixed with the current day's NPP share available for BNF.

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{N}_{\text{fix}}& {\displaystyle }\left(\text{NPP}\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{N}_{\text{fix,need}},& \\ & \mathrm{\text{if }}{\mathrm{\text{cost}}}_{\text{BNF}}\cdot N_{\text{fix,need}}<f_{\text{fixer}}\cdot f_{\text{NPP}}\cdot \text{NPP}\\ f_{\text{fixer}}\cdot f_{\text{NPP}}\cdot \text{NPP}/{\text{cost}}_{\text{BNF}},& \mathrm{\text{otherwise}},\end{array}\right.\end{array}\end{array}$$

where f_NPP is the maximum share (dimensionless) of NPP available for BNF, which is set to 0.14 (Kull, 2002) for the natural PFTs and to 0.25 for soybean and pulses. The average N fixer fraction (f_fixer) is set to 0.05 for the tropical, to 0.01 for the temperate and to 0.03 for the boreal zone (Yu and Zhuang, 2020). PFTs only fix additional N if the N uptake from other sources is insufficient and the net primary productivity (NPP) is larger than zero. The costs of BNF are set at a moderate constant value of 6 g C g N⁻¹ (Boote et al., 2009; Ryle et al., 1979; Patterson and Larue, 1983; Kaschuk et al., 2009).

2.4 Modelling protocol

To compare the two BNF approaches, we simulated two scenarios. The first is a potential natural vegetation (PNV) scenario, which does not include anthropogenic land use or agricultural production systems. The second is a scenario that includes agricultural land use (LU). The same input data sets were used for all scenarios. We used the climate data from the GSWP3-W5E5 data set (Kim, 2017; Cucchi et al., 2020; Lange et al., 2022), historical atmospheric N deposition (Yang and Tian, 2020), historical atmospheric CO₂ concentrations (Büchner and Reyer, 2022), historical land-use patterns (Ostberg et al., 2023) and grazing management data (Stenzel et al., 2024). For both BNF approaches, we conducted spin-up simulations of 3500 years using a random permutation of the climate data from 1901 to 1930. These spin-up simulations ensure that the C and N balances are in equilibrium. Afterwards, land use is introduced, and a second spin-up period of 390 years is run to capture the effects of historical land-use change on the C and N cycle. Following the two spin-up simulations, the model is run from 1901 until 2016 using the transient input data.

2.5 Model evaluation

We compared simulated total global BNF for both approaches against several estimates which were derived empirically or reported in other modelling studies. Data on these estimates are available from Davies-Barnard and Friedlingstein (2020a). The global BNF is calculated as the sum of BNF per area times grid cell area over all grid cells:

$$\begin{array}{}\text{(8)}& {\text{BNF}}_{\text{glob}}=\sum _{\text{cell}}^{n_{\text{cell}}}{\text{BNF}}_{\text{cell}}\cdot {\text{area}}_{\text{cell}}.\end{array}$$

For the evaluation we calculate the median, minimum and maximum between 2001 and 2010 and qualitatively compare these values against past estimates. We calculated the overlap between our results and the reported data if minimum and maximum values were available.

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }\text{Overlap}=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{0}& \mathrm{\text{if }}{x}_{\text{min}}>y_{\text{max}}\text{ or }{x}_{\text{max}}<y_{\text{min}}\\ (min(x_{\text{max}},y_{\text{max}})-max(x_{\text{min}},y_{\text{min}}\left)\right)/(y_{\text{max}}-y_{\text{min}})& \\ & \mathrm{\text{otherwise}},\end{array}\right.\end{array}\end{array}$$

where x_min and x_max are the simulated minimum and maximum, and y_min and y_max are the minimum and maximum values from the literature.

In addition, we compared our results to data obtained at several sites for the natural vegetation (Davies-Barnard and Friedlingstein, 2020a) and legume crops (Ma et al., 2022). To evaluate legume crop BNF and yields, we conducted additional local simulations matching the coordinates of the experiments following the protocol described in Sect. 2.4 but ensured that the respective crops (soybean or pulses) were grown under the reported water management (rainfed or irrigated). We calculated the root mean square error (RMSE) as follows:

$$\begin{array}{}\text{(10)}& \text{RMSE}=\sqrt{\sum _n^N(x_n-y_n{)}^{\mathrm{2}}/N},\end{array}$$

where N is the number of observations, and x_n and y_n are the simulated and observed values.

3 Results

3.1 Comparison of the BNF approaches

Comparing the simulated BNF of both approaches to data from the literature and experiments showed a substantial improvement of the global BNF (Sect. 3.1.1), as well as the latitudinal and spatial patterns (Sect. 3.1.2).

3.1.1 Comparison to data and other models

The two approaches show large differences in the simulated BNF. While the median global BNF between 2001 and 2010 was 191 Tg N yr⁻¹ for the Original approach, for the C-costly approach it was substantially lower, with a value of 109 Tg N yr⁻¹ (Fig. 2a). Comparing the global BNF of both approaches to estimates from the scientific literature shows the agreement of the C-costly values with several other data sources, while the Original approach overestimates most of the literature values. In particular, the recent estimate by Davies-Barnard and Friedlingstein (2020a) was closely matched by the C-costly approach, and 60 % of the simulated data were within the range of the Davies-Barnard and Friedlingstein (2020a) data (Fig. 2a). Despite the fact that the Original approach was not derived from the empirical relationship of Cleveland et al. (1999) for the legume crops, the data from Cleveland et al. (1999) are well matched by the Original approach, and only the spread of the Cleveland et al. (1999) data is underestimated. In comparison to the data of Xu-Ri and Prentice (2017), who reported much higher values compared to the other studies, BNF is underestimated by both approaches implemented in LPJmL. However, large differences are to be expected, considering that their approach does not calculate the actual BNF but rather the BNF needed to sustain global NPP (Xu-Ri and Prentice, 2017).

Figure 2 Evaluation against global (a) and site-specific data (b, c). Global evaluation plot inspired by Davies-Barnard and Friedlingstein (2020a) showing global BNF (in Tg N yr⁻¹) from different studies (black) compared to the Original (red) and C-costly (blue) BNF approach implemented in LPJmL. Studies are labelled by author names and whether they consider potential natural vegetation (PNV), actual natural vegetation (NV) or actual land use (LU). We assigned the Davies-Barnard and Friedlingstein (2020a) data to the LU category because they consider cropland area as grasslands and not as potential forest areas. Percentage values give the overlap between the ranges of the simulation results and the literature estimates derived using Eq. (9). Simulated values are the median between 2001 and 2010, and the ranges show the minimum and maximum. Site-specific evaluation (b, c) comparing data from observations for soybean (green) and pulses (blue) for rainfed (RF) (circle) and irrigated (IR) (triangle) experiment and simulation results is shown using the Original (b) and C-costly (c) BNF approaches. Labels show the RMSE of the two approaches.

Comparing the spatial patterns of the two approaches to those of Davies-Barnard and Friedlingstein (2020a) shows that the Original approach generally overestimated BNF in large areas of the tropics and temperate zones (Fig. C4c). The C-costly approach still overestimates BNF in the tropics and the production areas of soybean and/or pulses in India and the United States of America (USA), but the values are substantially smaller than in the Original approach (Fig. C4f). In both approaches, observed BNF is slightly underestimated in the central to western part of the USA, Canada, China, Kazakhstan, Russia and Mongolia.

On croplands, BNF was 21 Tg N yr⁻¹ with the C-costly approach, which is within the range of 17 to 31 Tg N yr⁻¹ reported by a recent review (Zhang et al., 2021) and other studies (Bodirsky et al., 2012; Chang et al., 2021). This contrasts with the overestimation of cropland BNF in the Original approach, which was 68 Tg N yr⁻¹. For the two legume crop functional types representing soybean and pulses, we compared the simulation results to BNF and yield data from experiments (Figs. 2b and c and C1a and b). For all except two experiments, the Original approach strongly overestimated BNF independent of the crop and the irrigation management. Using the C-costly approach, the cropland BNF was strongly reduced by a factor of approximately 2, leading to substantially lower RMSEs. While simulation results were closer to observations, some deviations remain. Pulses generally showed lower BNF for both approaches compared to soybean, while irrigated simulations generally showed a higher BNF and overestimated BNF compared to observations for all experiments in the Original approach and for the vast majority in the C-costly approach. Crop yields barely differed between the two approaches and were comparable to observations (Fig. C1a and b).

3.1.2 Global variation in BNF

Generally, BNF decreases from low to high latitudes with similar gradients but from different levels for the two approaches (Fig. 3). In latitudes with a high share of crop legumes (e.g. 30 to 40° S), the reduction in the BNF in the C-costly approach is especially large. Both the Original and the C-costly approaches underestimate BNF at high latitudes (the Original more strongly so) compared to Davies-Barnard and Friedlingstein (2020a). The C-costly approach shows a good performance in the mid-latitudes, but both approaches overestimate BNF compared to observations in the tropics (Fig. 3). In the Original approach, especially the high BNF of cropland contributes to the overestimation. For the low latitudes, both approaches exceed the values from Davies-Barnard and Friedlingstein (2020b). However, the higher BNF in the tropics is comparable to the median of the TRENDY-N ensemble (see Sect. 4 and Kou-Giesbrecht et al., 2023).

Figure 3 Simulated average annual BNF (in $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) for the years 2001 to 2010 using the Original (a) and C-costly (b) approaches. The average (line) and 5th to 95th percentiles (shading) of simulated and observed BNF per latitude (in $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), using the Original (red) and C-costly (blue) approaches and data from Davies-Barnard and Friedlingstein (2020a) (DBF), are shown (c).

With the Original approach, mineral N was added to the first soil layer and subsequently incorporated by the PFTs via the passive and active N-uptake pathway. This did not allow a separate identification of N to be taken up via BNF from the total N uptake, except for the legume crops which fixed their entire N deficit. Using the C-costly BNF, the model separates N uptake by BNF from passive and active N uptake against N concentration gradients (Marschner et al., 1991; Fisher et al., 2010), facilitating the analysis of the share of BNF in total N uptake subsequently referred to as %N_dfa, which is commonly used to refer to this variable in the empirical literature (e.g. Herridge et al., 2008). In the PNV simulation, values for %N_dfa were between 0 % and 20 % for most of the grid cells (Figs. C7 and C8b). The distribution is bimodal, showing a peak below 5 % and one at approximately 10 %. For the dynamic land-use simulation, the values for %N_dfa are similar, but the second peak is barely distinguishable because of a higher share of %N_dfa values between 5 % and 10 % (Figs. C7 and C8a). For the crop legumes, %N_dfa is substantially higher, with the peak at around 80 % (Fig. C7d). The highest values are simulated at low latitudes, especially in India, sub-Saharan Africa and South America, while the lowest values are simulated in Canada, Russia and southern China (Fig. C8d). In the Original approach, %N_dfa of legume crops was almost exclusively 100 % (Fig. C7e).

In both approaches, BNF per area is higher for agricultural land than for natural vegetation (Fig. C3d and f). BNF is especially high in hot spots of legume crop production such as Argentina, Brazil, India and the USA (Fig. 3a and b). While the spatial pattern is similar between the two approaches, in the Original approach, the cropland BNF leads to prominent peaks in the latitudinal distribution (Fig. 3c). These peaks correspond to hot spots of legume crop production, where the C-costly approach is up to 15 $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ lower (Fig. C4).

For natural vegetation, the differences are smaller. Here, the BNF in the Original approach is up to 4 $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ higher compared to the C-costly approach (Fig. C4). Here, the spatial patterns differ and show a stronger reduction in the BNF in dry regions (e.g. central Australia, the Eurasian steppe regions, southeast China and parts of Africa).

The various natural PFTs contribute differently to the lower overall BNF in the C-costly approach (Figs. C5 and C6). To some extent, this reflects changes in the PFT distribution (Figs. S1 and S2 in the Supplement). For the tropical PFTs, BNF is lower for the broadleaved raingreen tree (Δ5.25 Tg N yr⁻¹; Fig. C5b) and the herbaceous PFT (Δ14.1 Tg N yr⁻¹; Fig. C5i) and higher for the broadleaved evergreen tree (Δ7.3 Tg N yr⁻¹; Fig. C5a). While the temperate needleleaved evergreen tree PFT contributed to BNF in low latitudes outside its expected habitat (e.g. in India and Brazil) in the Original approach, this issue was resolved with the C-costly approach (Fig. C5c). The temperate PFTs all fix less N in the C-costly approach than in the Original approach. The reductions are smaller for the broadleaved evergreen (Δ3.6 Tg N yr⁻¹; Fig. C5d) and summergreen (Δ3.8 Tg N yr⁻¹; Fig. C5e) tree and the herbaceous PFT (Δ4.7 Tg N yr⁻¹; Fig. C5j) compared to the needleleaved evergreen tree (Δ9.1 Tg N yr⁻¹; Fig. C5c). The BNF of boreal PFTs is similar (Δ around 0.5 Tg N yr⁻¹; Fig. C5f, g, and k) for all PFTs, except the needleleaved summergreen tree (Δ1.2 Tg N yr⁻¹; Fig. C5h), which fixes less N in the C-costly approach. In the Original approach, the temperate herbaceous PFT contributed twice as much as in the C-costly approach to the biological N fixation of the polar vegetation (Fig. C5j). For pulses, the BNF was 14.6 Tg N yr⁻¹, and for soybean,the BNF was 6.4 Tg N yr⁻¹ lower with the C-costly approach.

3.2 Effects on the nitrogen and carbon cycle and productivity

In LPJmL the C and N cycles are coupled via, for example, the N limitation of gross primary productivity (GPP), which controls the amount of assimilated C, the role of plant organ carbon-to-nitrogen (C:N) ratios for maintenance respiration and the availability of the resulting NPP for BNF. Additionally, the N content of the different plant organs (leaves, roots, sapwood, heartwood and storage organs) is derived dependent on the respective C content ensuring that their C:N ratios remain within a prescribed range. As a result, the N balance components presented in the following section are strongly shaped by their C cycle counterparts as the overall C and N balances represented by LPJmL are intimately linked.

We describe the N balance as the sum over in- and outfluxes of the vegetation and the soil. Therefore, the overall balance contains a change in vegetation and soil N stocks, including organic and mineral forms of N.

3.2.1 Potential natural vegetation

Simulating only natural vegetation resulted in a positive terrestrial N balance with an average sink of 52 Tg N yr⁻¹ for the Original approach and 54 Tg N yr⁻¹ for the C-costly approach between 2001 and 2010 (Fig. 4a and b and Table 2). In 1901, N in- and outputs were almost balanced, and the sink remained small until the 1950s when N inputs from deposition increased, resulting in an increased sink. While the overall N balance was similar for both BNF approaches, the size of several components was different. The total BNF simulated with the Original approach was approximately double that of the C-costly BNF, leading to higher soil mineral N and organic C and N stocks. However, mineral N stocks were not utilised by the vegetation but instead lost to the atmosphere and waterbodies, leading to higher N emissions and leaching using the Original approach. Here, 112 Tg N yr⁻¹ were emitted and 56 Tg N yr⁻¹ were leached on average between 2001 and 2010, while for the C-costly approach, only 79 Tg N yr⁻¹ were emitted and 39 Tg N yr⁻¹ were leached (Table 2). All types of emissions are lower with the C-costly approach. NH₃ emissions from volatilisation decrease by 14 Tg N yr⁻¹, N₂ emissions by 12 Tg N yr⁻¹, N₂O emissions by 3 Tg N yr⁻¹ and fire emissions by 5 Tg N yr⁻¹. Synchronised with the increase in deposition over time, emissions and leaching also increase in both approaches, with stronger increases in the C-costly approach. Overall, N inputs increased by 35 Tg N in total from 1950 to 2000 in the Original approach and by 42 Tg N in the C-costly approach, while N losses from emissions and leaching increased by 1 and 4 Tg N, respectively (Table 2).

Figure 4 Global terrestrial N balance. Scenarios include the Original approach, the C-costly approach for natural vegetation and actual land use. Net balance is denoted by the black line. N inputs include N from manure, synthetic fertiliser deposition, PFT establishment (Estab) and BNF. N losses include leaching, volatilisation, N₂ emissions, fire N, harvested N, land-use change emissions (deforestation and product turnover) and N₂O emissions from nitrification and denitrification.

Table 2 N balance values for 2001 to 2010 shown in the figures. LUC (land-use change) includes deforestation emissions, product turnover and negative N fluxes.

^a Bouwman et al. (2013). ^b Zaehle et al. (2010). ^c Braakhekke et al. (2017). ^d Volatilisation from natural soils (2.4 Tg N yr⁻¹) from Bouwman et al. (2002). ^e Volatilisation from manure and synthetic fertiliser on croplands and grasslands (19 Tg N yr⁻¹) from Bouwman et al. (1997). ^f Galloway et al. (2004). ^g Sutton et al. (2013). ^h Tian et al. (2019). ⁱ Tian et al. (2020). ^j Scheer et al. (2020). ^k Yu and Zhuang (2020).

In addition to the changes in several N cycle components, we excepted changes in C cycle components. Overall, the C balance was similar for both approaches (Fig. C9a and b). For the PNV simulations, the only C input into the system was the NPP. The NPP was 2.2 Pg C yr⁻¹ lower with the C-costly approach compared to the Original approach. However, C losses from heterotrophic respiration and fire were also lower by 1.9 and 0.3 Pg C yr⁻¹, respectively.

3.2.2 Dynamic land use

The simulations with dynamic land use include agricultural production and related additional N in- and outputs. Additional inputs are N from the application of manure and synthetic fertilisers, and additional outputs are N removed through crop harvesting, grazing and emissions from land-use change. The differences in the total BNF, soil mineral N and organic C and N stocks are similar to the PNV simulations. Between 2001 and 2010, LPJmL simulated an average N sink of 44 Tg N yr⁻¹ for the Original and 45 Tg N yr⁻¹ for the C-costly approach (Fig. 4c and d and Table 2). Already in 1901, the N balances of the PNV and dynamic land-use simulations diverge. Since there are no synthetic fertiliser inputs in 1901, only the relatively small additional inputs from establishment and manure were counteracted by N removal through crop harvesting and land use change emissions, shifting the total N balance towards a smaller source. This persists even after inputs from manure and fertiliser were increased from the 1950s onwards, which not only resulted in higher crop yields, and therefore N removed through harvesting, but also increased N losses from emissions and leaching. As was shown for the PNV simulations, the overall N balance is similar for both approaches but with different in- and output terms driven by the higher BNF in the Original approach. N emissions and leaching for the Original approach (128 and 74 Tg N yr⁻¹, respectively) were higher than for the C-costly approach (99 and 55 Tg N yr⁻¹, respectively) (Table 2). Similar to the PNV scenario, all types of emissions are lower with the C-costly approach. NH₃ emissions from volatilisation decrease by 11 Tg N yr⁻¹, N₂ emissions by 13 Tg N yr⁻¹, N₂O emissions by 3 Tg N yr⁻¹ and fire emissions by 2 Tg N yr⁻¹. N removal from harvesting was 142 Tg N yr⁻¹ on average between 2001 and 2010 for the Original approach and 108 Tg N yr⁻¹ for the C-costly approach. This indicates a stronger N limitation of agricultural areas in the C-costly approach. The majority of this reduction can be attributed to managed grassland and not croplands (Figs. S3 and S4).

Similar to the PNV simulations, the overall C balance barely differed between the two approaches (Fig. C9c and d). While the C input from manure and establishment was similar for both approaches, NPP was 1.8 Pg C yr⁻¹ lower in the C-costly approach. C lost from land-use change was similar. Fire emissions and C removed through harvesting only differed by 0.1 Pg C yr⁻¹, while heterotrophic respiration was 1.4 Pg C yr⁻¹ lower in the C-costly approach than in the Original approach.

4 Discussion

While the Original approach only indirectly accounts for temperature and water limitation, as these also limit evapotranspiration and NPP, the C-costly approach explicitly considers the limitation of BNF from soil temperature, water content and NPP separately. These have long been established as limiting factors for BNF. Depending on the prevailing conditions, BNF may be limited more strongly by temperature or soil moisture or a combination of both. The role of temperature was explored early on by Meyer and Anderson (1959) and was followed by numerous studies for different plant species or legume crop varieties and temperature ranges (e.g. Montañez et al., 1995). Such studies allowed the quantification of optimal temperature ranges and limits facilitating the development of functions such as f_T(T_soil), which are based on empirically derived temperature curves (e.g. Halliday and Pate, 1976). A similarly large literature body exists on the role of soil moisture (e.g. Serraj et al., 1999; Rousk et al., 2018). Valentine et al. (2018), describing the different pathways with respect to how drought stress inhibits BNF, with an important aspect being the change in nodule water potential that indicates a strong connection to soil water content. While flooding of soils can also inhibit BNF through O₂ limitation, nitrogenase was shown to be more active in waterlogged environments (Jiang et al., 2021). Therefore, we are confident that our linear function for f_W(SWC), assuming that only soil moisture levels that are too low limit BNF (McGechan, 1999), reflects empirical observations well. As BNF is associated with a respiratory loss of C, the net amount of C assimilated via photosynthesis (NPP) available for BNF, as well as the fixation efficiency (respiratory loss of C per gained N), forms an additional important controlling factor. A recent meta-analysis by Yao et al. (2024) highlights the importance of plant taxa for BNF in addition to abiotic factors. This is in line with early experimental work that quantified the respiratory loss of C per N fixed (Reed et al., 2011; Patterson and Larue, 1983; Voisin et al., 2003) and the total amount of NPP spent on fixation (Kaschuk et al., 2009) for different N-fixing plants and already showed that functional traits have to be considered when assessing BNF. Therefore, including NPP and a cost of fixation, as we did with the C-costly approach, is an important conceptual improvement.

The C-costly approach is not only conceptually superior to the simplistic Original approach in LPJmL, it also performs better in comparison to external data. Still, some mismatches with reference data remain, such as an overestimation of BNF in the tropics (Fig. 3c). However, the ensemble mean of a recent study evaluating the N cycle of 11 DGVMs shows a similar overestimation in the tropics and a large bias, indicating little agreement between models (Kou-Giesbrecht et al., 2023). They attributed this to the fact that BNF is typically modelled as a function of vegetation activity expressed either through NPP or evapotranspiration. Our results show that the overestimation of tropical BNF is reduced if temperature and water availability are considered separate limitations, which supports their interpretation. Furthermore, the NPP that can be used for BNF depends on the overall productivity, which certainly is higher in the tropics. It is likely that additional variables not considered in our approach constrain BNF there, such as phosphorus limitation (Vitousek, 1984; Lee et al., 2019). However, it has also been suggested that as a result of higher N losses, tropical BNF should be higher than observations imply (Hedin et al., 2009). This could be a result of uncertainties inherent to BNF measurements (Soper et al., 2021) or the limited amount of data available from tropical ecosystems.

Furthermore, simulated BNF was at the higher end of the range reported by Davies-Barnard and Friedlingstein (2020b) for the C-costly approach. One explanation is that Davies-Barnard and Friedlingstein (2020b) aggregate cropland and grassland areas, assuming that their BNF rates are identical. However, a recent study provides evidence that the BNF of crop legumes might actually be substantially higher than that of forage legumes (Herridge et al., 2022; Peoples et al., 2021), and therefore, the BNF of croplands and grassland cannot be assumed to be similar. Consistent with this, we also had to select much higher potential N fixation rates for the crop PFTs compared to the other PFTs to achieve sufficient cropland BNF (Table 1).

We expected that limiting BNF of legume crops would result in stronger N stress and reduced yields. However, yields for the legume crops were similar between the two approaches. One explanation is the direct link that the maintenance respiration of a plant organ has to its N content. Reducing the N that is taken up via BNF results in a lower organ N content and maintenance respiration and thus similar NPP. Indeed, C:N ratios are higher for the C-costly approach compared to the Original approach, indicating a lower plant N content (Fig. C2).

The average %N_dfa of legume crops was between approximately 30 % and 100 % for the C-costly approach and 100 % for the Original approach. Comparing the distribution (Fig. C7d) to %N_dfa observations shows that the values of the C-costly approach are possible but at the upper end of observations, while those of the Original approach are not supported by observations. For soybean, experimental values range from 0 % to 98 %, with an average of 52 % (Salvagiotti et al., 2008). Herridge et al. (2008) report average values between 40 % and 75 % on average and up to 97 % for experiments but only 36 % and 68 % for farmers' fields, depending on the cultivated legume crop. %N_dfa is strongly related to soil mineral N content and thus fertilisation levels. The high %N_dfa may be an indication that either fertiliser levels or active and passive N uptake and retranslocation of N at leaf senescence are underestimated by LPJmL, and respective processes should be re-evaluated. We found a higher %N_dfa for both the natural vegetation and the cropland in warm and dry areas (Fig. C8) where mineralisation of organic N is limited (Dessureault-Rompré et al., 2010).

We expected that the differences in the BNF between the two approaches would be reflected by differences in the C stocks and fluxes due to the close link of the C and N cycles in LPJmL. Both the C inflow into terrestrial C stocks from NPP and outflows from harvest, heterotrophic respiration and fire were lower in the C-costly approach, leading to a similar net C balance for the two approaches (Figs. C9 and S5). Accounting for the cost of BNF in the form of respiratory losses of NPP leads to lower NPP, which limited biomass accumulation and in turn harvest, as well as biomass available for burning and heterotrophic respiration via reduced litter accumulation. Because of the close link of the C and N cycles, the net N balance is also similar for the two approaches. The lower BNF in the C-costly approach results in lower N outfluxes, i.e. leaching, emissions and harvests. The Original approach added mineral N to the soils of the natural vegetation and managed grassland even if the vegetation was not N limited. Legume crops that received all the N that they demanded, as in the Original approach, returned high N content residues to the soil, increasing N inputs and mineral N stocks. As a result, the mineral N content of soils was higher in the Original approach, explaining the differences in yields and leaching. Similarly, soil mineral N content influences N emissions except fire emissions, which are controlled by the N content of the burned vegetation and litter. Since this also decreased, fire emissions were lower with the C-costly approach. In contrast to the lower BNF, which is in line with observations, N losses from leaching and emissions (from volatilisation, denitrification, nitrification, fire and land-use change) are underestimated by LPJmL simulations compared to observational data (see Table 2) in both approaches. The overestimation of emissions from volatilisation of soil ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ is strongly reduced with the C-costly approach because the soil ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool is lower in the C-costly approach compared to the Original approach, where BNF is directly added to the soil ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ (see Sect. 2.2). While N₂O emissions compare well to literature estimates, N₂ emissions are more strongly underestimated with the C-costly approach. Similar to the soil ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool, the soil ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ is reduced because less ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ is available for nitrification, resulting in reduced N₂ emissions. Overall, the reduction shifts emissions from an over- to an underestimation of the literature values. While one source of differences is the missing representation of NO_x emissions in LPJmL, this is not sufficient to fully explain the difference. However, the models of the TRENDY-N ensemble also underestimated N losses from emissions of NH₃, N₂O, NO_x, and N₂, as well as leaching (Kou-Giesbrecht et al., 2023), suggesting that processes within DGVMs and scenario assumptions need to be revised. For LPJmL, we identified several potential causes. First, the manure input accounts only for manure applied to cropland, and the total amount is in line with other sources reporting cropland manure (Zhang et al., 2021), but does not account for manure added to grasslands other than the internal recycling by grazing animals (Heinke et al., 2023). Second, N losses and emissions strongly vary between different agricultural production systems whose representation would require not only the implementation of more detailed management options but also data sets on the spatial patterns of the application of different management specifics of these systems. Third, we conducted our simulations assuming cover cropping outside the growing season on all croplands, which overestimates the extent of cover cropping and reduces N losses. However, data on cover cropping systems are not available (e.g. Porwollik et al., 2022).

While the C-costly approach improved simulation results for BNF, as well as other components of the N balance, and model results are in line with other DGVMs that represent the N cycle, we see potential for further improvement. The C-costly approach depends on multiple parameters, some of which are not well constrained. Values for the potential N fixation rate vary between species and across sites (Ma et al., 2022), and selecting one value to be representative of one PFT or even all PFTs of an entire climate zone is a strong simplification. Furthermore, experiments have shown a large variation in the respiratory cost of BNF (Reed et al., 2011; Patterson and Larue, 1983; Voisin et al., 2003), as well as the amount of NPP that different plant species invest (Kaschuk et al., 2009), which is not well reflected by the current parameterisation.

In addition, we assume a constant fraction of N fixers present in a community. However, the number of N fixers changes over time, dependent on N stress (Herben et al., 2017; Taylor et al., 2019). N fixation, the share of fixers and/or nodule abundance is low in undisturbed N-rich environments, and nodules need to be produced to increase N fixation if N availability decreases (Fisher et al., 2010; Crews, 1999). Similarly, N fixation does not cease instantaneously when N becomes more abundant but is only reduced after the share of fixers and/or nodule abundance has decreased (Thornley et al., 1995; Herben et al., 2017). In contrast, fixers are always present in LPJmL and can instantly fix N if necessary. Therefore, LPJmL likely simulates an adaptation response that is too quick for the changing N availability and overestimates the short-term capability of the community to buffer changes in N availability.

While our approach simulates the total amount of BNF well, it does not distinguish the symbiotic from free-living or heterotrophic N fixation. However, these are two different sources of N, and their share of total BNF shows large spatial heterogeneity (Davies-Barnard and Friedlingstein, 2020b). In contrast to symbiotic BNF, free-living BNF does not require NPP expenditures, and separating the two may further improve simulation results for NPP and its dependent variables.

In the following, we qualitatively compare our approach to common approaches used in crop models and DGVMs. A synthesis of nine crop models by Liu et al. (2011) showed that soil water status and N supply were the most widely considered control variables. Soil temperature was only considered by four models and plant C supply only by two models, despite their importance for limiting BNF. The C-costly approach also uses empirical factors to account for soil temperature and soil water status, whereas the role of N supply, plant C supply and plant growth stage are simulated mechanistically in LPJmL, which is a clear distinction from the models assessed by Liu et al. (2011).

Our approach is at the higher end of the complexity when compared to 11 TRENDY-N DGVMs that include the N cycle. As shown in Kou-Giesbrecht et al. (2023), five DGVMs follow an approach similar to the Original approach, calculating BNF based on evapotranspiration or NPP; three models calculate BNF as a function of N limitation; two models assume a constant BNF; and in one model BNF is derived in post-processing to close the N cycle. The remaining three models use more complex approaches which can be compared to the C-costly approach. The Community Land Model (CLM) version 5.0 (Lawrence et al., 2019) uses an approach based on Fisher et al. (2010), explicitly minimising the cost of active N uptake, retranslocation and BNF and distinguishes asymbiotic from symbiotic N fixation. In the Canadian Terrestrial Ecosystem Model (CTEM), BNF is a function of temperature, vegetation cover, soil nitrate and plant structural C pools (Arain et al., 2006; Dickinson et al., 2002). The Dynamic Land Ecosystem Model (DLEM) considers soil temperature, soil moisture status, soil C and soil N (Tian et al., 2015). While the approach used in CLM5.0 is more complex compared to the C-costly approach and addresses some of the conceptual shortcomings of the C-costly approach discussed earlier, the approach used in CTEM is of similar complexity and simulates values at the upper end of the recent literature estimates (Kou-Giesbrecht et al., 2023). However, global BNF values and latitudinal distribution simulated by CLM5.0, as shown by Kou-Giesbrecht et al. (2023) in Figs. 3 and A6, are comparable to those simulated with C-costly approach. To fully assess the advantages of such a complex approach over the C-costly approach or that of CTEM or DLEM, a comparison of the spatial patterns or of simulations at higher spatial resolution could be a worthwhile future endeavour.

5 Conclusions

Compared to the simplistic Original BNF implementation in LPJmL, the more complex C-costly approach, as described here, presents a substantial improvement of the representation of BNF in LPJmL. While the Original approach led to an overestimation of BNF and was insensitive to soil temperature and soil water conditions, the C-costly approach overcomes these issues and can help to better project future BNF and its effects on N limitation of the terrestrial biosphere, as well as losses of reactive N to the environment, including the greenhouse gas nitrous oxide (N₂O). Further research is needed, especially with respect to balancing different in- and outfluxes and internal recycling rates. The current improvement of BNF simulations with LPJmL and the associated underestimation of loss terms exemplifies the scope of this problem. Our study highlights the importance of a detailed implementation of the processes controlling BNF for N cycling in DGVMs. While the C-costly approach already improved simulation results, we think that additional benefits would be gained by explicitly separating BNF by symbiotic and free-living bacteria and from accounting for the costs of other N-uptake sources, except passive N uptake.

Appendix A: Nitrogen demand and uptake

The total N demand (N_demand in g N m⁻²) at any time t is the sum of the leaf N demand for RuBisCo and structural components (N_demand,leaf in g N m⁻²) and the N demand for structural components of the other plant compartments.

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{N}_{\text{demand,leaf}}=& {\displaystyle }\mathrm{25}\cdot \mathrm{0.02314815}/\text{daylength}\cdot V_{\text{max}}\\ {\displaystyle }& {\displaystyle }\cdot \exp (-\mathrm{0.02}\cdot (T-\mathrm{25}\left)\right)\cdot f_{\text{LAI}}\left(\text{LAI}\right)\\ {\displaystyle }& {\displaystyle }+N{C}_{\text{leaf,median}}\cdot C_{\text{leaf},t},\end{array}\end{array}$$

where V_max (g C m⁻²) is the PFT-specific maximum carboxylation capacity computed based on absorbed photosynthetically active radiation (APAR) and canopy conductance (Schaphoff et al., 2018b; Sitch et al., 2003). T is the average temperature (°C) of the current day, and daylength is the duration of daylight (h). f_LAI(LAI) is a dimensionless modifier to account for the current leaf area index (von Bloh et al., 2018), and C_leaf (g C m⁻²) is the current leaf C content.

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{f}_{\text{LAI}}& {\displaystyle }\left(\text{LAI}\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}max(\mathrm{0.1},\text{LAI})& \text{for }\text{LAI}<\mathrm{1}\\ \exp (\mathrm{0.08}\cdot min(\text{LAI},\mathrm{7}\left)\right)& \text{otherwise,}\end{array}\right.\end{array}\text{(A3)}& {\displaystyle }{C}_{\text{leaf},t}=C_{\text{leaf}}+f_{\text{leaf}}\cdot \sum _{t^{\prime }=\mathrm{1}}^t{\text{NPP}}_{t^{\prime }}-\mathrm{\Delta }{\text{litter}}_{t^{\prime }}.\end{array}$$

LAI is the current leaf area index, and ${\sum }_{t^{\prime }=\mathrm{1}}^t{\text{NPP}}_{t^{\prime }}-\mathrm{\Delta }{\text{litter}}_{t^{\prime }}$ is the difference between the accumulated biomass increment and litterfall.

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{N}_{\text{demand},t}=\\ {\displaystyle }& {\displaystyle }(N_{\text{demand,leaf}}+{\sum }_{m=\mathrm{1}}^M{N}_m+{\text{NC}}_t\cdot {\sum }_{m=\mathrm{1}}^M(f_m/R_m)\\ {\displaystyle }& {\displaystyle }\cdot {\sum }_{t^{\prime }=\mathrm{1}}^t{\text{NPP}}_{t^{\prime }}-\mathrm{\Delta }{\text{litter}}_{t^{\prime }})\cdot (\mathrm{1}+k_{\text{store}}),\end{array}\end{array}$$

where M is one for grasses, two for trees and three for crops, equalling the number of the respective PFT plant compartments and excluding leaves, $N{C}_t=min(max(N_{\text{leaf},t}/C_{\text{leaf},t},N{C}_{\text{leaf,low}}),N{C}_{\text{leaf,high}})$. f_m is the fraction of biomass allocated to the compartment m, R_m is the C:N ratio of compartment m relative to the leaf C:N ratio, and k_store is a PFT-specific parameter to maintain the PFT labile N storage. Passive and active N uptake (N_uptake) from each soil layer l (n_soillayer=6) is calculated as a function of the potential N uptake of the root system.

$$\begin{array}{}\text{(A5)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{N}_{\text{uptake}}={\sum }_{l=\mathrm{1}}^{n_{\text{soillayer}}}\mathrm{2}& {\displaystyle }\cdot N_{\text{up,root}}\cdot C_{\text{root},t}\cdot {\text{rootdist}}_l\\ {\displaystyle }& {\displaystyle }\cdot f_N\left(N_{\text{avail},l}\right)\cdot f_T\left(T_{\text{soil},l}\right)\\ {\displaystyle }& {\displaystyle }\cdot f_{\text{NC}}\left({\text{NC}}_{\text{plant}}\right),\end{array}\end{array}$$

where N_up,root is the PFT-specific maximum N-uptake rate per unit of fine root mass in each layer, C_root,t is the current root C, rootdist_l is the fraction of roots in layer l. f_N, f_T and f_NC are dimensionless modifiers for the availability of mineral N, soil temperature and plant N:C ratio (von Bloh et al., 2018). C_root,t is calculated as C_leaf,t in Eq. (A3). The root distribution can be calculated from the proportion of roots from the surface to soil depth z, following Jackson et al. (1996):

$$\begin{array}{}\text{(A6)}& {\text{rootdist}}_z={\displaystyle \frac{\mathrm{1}-{\mathit{\beta }}_{\text{root}}^z}{\mathrm{1}-{\mathit{\beta }}_{\text{root}}^{z_{\text{bottom}}}}},\end{array}$$

where z_bottom is the lower boundary of the last soil layer and β_root is a PFT-specific parameter (Table A1). The root proportion of one soil layer can be calculated as

$$\begin{array}{}\text{(A7)}& {\text{rootdist}}_l={\text{rootdist}}_{z\left(l\right)}-{\text{rootdist}}_{z(l-\mathrm{1})}.\end{array}$$

Table A1 PFT-specific parameters used in N-demand and N-uptake calculations.

f_N follows Michaelis–Menten kinetics,

$$\begin{array}{}\text{(A8)}& f_N\left(N_{\text{avail},l}\right)=k_{N,\text{min}}+{\displaystyle \frac{N_{\text{avail},l}}{N_{\text{avail},l}+K_{N,\text{min}}\cdot {\mathrm{\Theta }}_{\text{max}}\cdot d_{\text{soil}}}},\end{array}$$

where $N_{\text{avail},l}={{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}}^{-}+{{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l}}^{+}$, k_N,min and K_N,min are the PFT-specific parameters describing the Michaelis–Menten kinetics. Θ_max is the soil-type-specific fractional pore space, and d_soil,l (dimensionless) is the soil layer depth (m).

f_T is the temperature function given by Thornley (1991)

$$\begin{array}{}\text{(A9)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{f}_T& {\displaystyle }\left(T_{\text{soil},l}\right)=\\ {\displaystyle }& {\displaystyle }max\left({\displaystyle \frac{(T_{\text{soil},l}-T_{\mathrm{0}})\cdot (\mathrm{2}\cdot T_m-T_{\mathrm{0}}-T_{\text{soil},l})}{(T_r-T_{\mathrm{0}})\cdot (\mathrm{2}\cdot T_m-T_{\mathrm{0}}-T_r)}},\mathrm{0}\right),\end{array}\end{array}$$

where $T_{\mathrm{0}}<T_r<\mathrm{2}\cdot T_m-T_{\mathrm{0}}$ has to be fulfilled. von Bloh et al. (2018) defined T_m=15 °C, T_r=15 °C and $T_{\mathrm{0}}=-\mathrm{25}$ °C, which leads to the maximum of one at temperatures of 15 °C and higher and non-zero values above −25 °C.

f_NC was taken from Zaehle and Friend (2010),

$$\begin{array}{}\text{(A10)}& f_{\text{NC}}=min\left(max\left({\displaystyle \frac{{\text{NC}}_{\text{plant}}-{\text{NC}}_{\text{leaf,high}}}{{\text{NC}}_{\text{leaf,low}}-{\text{NC}}_{\text{leaf,high}}}},\mathrm{0}\right),\mathrm{1}\right),\end{array}$$

where ${\text{NC}}_{\text{plant}}=\frac{N_{\text{leaf}}+N_{\text{root}}}{C_{\text{leaf}}+C_{\text{root}}}$, NC_leaf,min and NC_leaf,max are PFT-specific parameters extracted from the TRY database (Kattge et al., 2020) (Table A1).

The labile-N values N_labile,t are the current reserves which have accumulated via N uptake and retranslocation.

$$\begin{array}{}\text{(A11)}& N_{\text{labile},t}=N_{\text{labile}}+\sum _{t^{\prime }=\mathrm{1}}^t{N}_{\text{uptake},t^{\prime }}-N_{\text{resorb},t^{\prime }}\end{array}$$

Appendix B: Spin-up simulation carbon stocks

With constant forcing (i.e. stable pre-industrial atmospheric CO₂ concentration, atmospheric N deposition and climate), the global C stocks showed a residual trend of −0.0106 Pg C yr⁻¹ for the Original approach and −0.0121 Pg C yr⁻¹ for the C-costly approach. This is 8–10 times lower than the steady-state criterion of the 0.1 Pg C yr⁻¹ residual trend after spin-up, which is used by the Global Carbon Project to validate DGVMs for inclusion in their global C budget analysis (Friedlingstein et al., 2022). At the grid cell level, the vast majority of cells (94 % for the Original approach and 95 % for the C-costly approach) exhibited residual trends in total C stocks of less than ±1 g C m² yr⁻¹. The corresponding maps are shown in Fig. B1.

Figure B1 Residual trends in C stocks after the spin-up simulation averaged over 1000 years for the Original approach (a) and the C-costly (b) approach.

Appendix C: Additional figures and tables

Figure C1 Simulated and observed crop yields (a, b) for soybean (green) and pulses (blue) and BNF in natural vegetation (c, d).

Figure C2 Vegetation C:N ratio for the years 2001 to 2010 for rainfed (RF) and irrigated (IR) soybean (red) and pulses (blue) for the Original approach and C-costly approach.

Figure C3 The 2001 to 2010 average BNF (in $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) of the potential natural vegetation (PNV) simulations (a, b), of the natural vegetation (NV) (c, d) and of the managed land (AG) (e, f) area fractions of the dynamic land-use (LU) simulations.

Figure C4 Difference between 2001 to 2010 average BNF (in $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) between the two approaches (a–d) for the dynamic land-use (LU) simulations (a), for the potential natural vegetation (PNV) simulations (b), for the area fractions of natural vegetation (NV) (c) and for the managed land (AG) (e) of the dynamic land-use simulations and the difference compared to the data from Davies-Barnard and Friedlingstein (2020a) (DBF) (e, f).

Figure C5 Latitudinal distribution of BNF for each PFT for the dynamic land-use simulations for the Original approach (red) and C-costly approach (blue).

Figure C6 Latitudinal distribution of BNF for each PFT for the potential natural vegetation simulations for the Original approach (red) and C-costly approach (blue).

Figure C7 Density distribution of the fraction of BNF of the total N uptake for the dynamic land-use simulations (a), for the potential natural vegetation (b) and for the area fractions of natural vegetation (NV) (c) and cropland (CL), using the C-costly approach (d) and the Original (e) approach for the dynamic land-use simulations.

Figure C8 Global distribution of the fraction of %N_dfa for the dynamic land-use (a), for the potential natural vegetation simulations (b) and the natural vegetation (c) and cropland (d) fraction of the dynamic land-use simulation.

Figure C9 Global terrestrial C balance. Scenarios include the Original approach and the C-costly approach for natural vegetation and actual land use. The net balance is denoted by the black line. C inputs include C from manure, PFT establishment (Estab) and NPP. C losses include heterotrophic respiration, fire emissions, harvested C and land-use change emissions (from deforestation and product turnover).