Implementing the nitrogen cycle into the dynamic global vegetation, hydrology, and crop growth model LPJmL (version 5.0)

The well-established dynamical global vegetation, hydrology, and crop growth model LPJmL is extended with a terrestrial nitrogen cycle to account for nutrient limitations. In particular, processes of soil nitrogen dynamics, plant uptake, nitrogen allocation, response of photosynthesis and maintenance respiration to varying nitrogen concentrations in plant organs, and agricultural nitrogen management are included in the model. All new model features are described in full detail and the results of a global simulation of the historic past (1901–2009) are presented for evaluation of the model performance. We find that the implementation of nitrogen limitation significantly improves the simulation of global patterns of crop productivity. Regional differences in crop productivity, which had to be calibrated via a scaling of the maximum leaf area index, can now largely be reproduced by the model, except for regions where fertilizer inputs and climate conditions are not the yield-limiting factors. Furthermore, it can be shown that land use has a strong influence on nitrogen losses, increasing leaching by 93 %.

1 Introduction

The dynamics of the terrestrial biosphere and the associated terrestrial carbon cycle are of central importance for Earth system science. Climate–carbon cycle feedbacks have become integral parts of Earth system models (ESMs) for climate change projections. However, the terrestrial carbon cycle dynamic are not only driven by climate and carbon dioxide (CO₂) fertilization (Schimel et al., 2015; Norby et al., 2005), but also by land use change (Müller et al., 2006, 2016; Arneth et al., 2017; Le Quéré et al., 2016) and vegetation dynamics (Müller et al., 2016, and references therein). Nutrient limitations, especially from nitrogen, are also important constraints on vegetation growth and the terrestrial carbon cycle. Smith et al. (2016) and Zaehle et al. (2015) suggested that Earth system models contributing to the CMIP5 data archive overestimate the response of net primary productivity to elevated CO₂ because the models largely miss the constraints from nutrient limitation. Also Wieder et al. (2015) find that nitrogen limitation may substantially reduce projected increases in net primary productivity (NPP) under climate change and elevated atmospheric CO₂ concentrations ([CO₂]), possibly even converting the terrestrial biosphere into a net carbon source by the end of the 21st century. Over the last decade, nitrogen limitation has been increasingly accounted for in dynamic global vegetation models (DGVMs) and ESMs (Thornton et al., 2007; Gerber et al., 2010; Zaehle et al., 2010b; Smith et al., 2014). The Lund Potsdam Jena managed Land (LPJmL) dynamic global vegetation, hydrology, and crop growth model has been widely applied to research questions on the terrestrial carbon cycle, hydrology, and agricultural production (Schaphoff et al., 2018b, and references therein) and performed similarly to other dynamic vegetation models (Friend et al., 2014; Warszawski et al., 2013; Chang et al., 2017), hydrology models (Schewe et al., 2014), and crop models (Müller et al., 2017). However, LPJmL so far has not explicitly accounted for nutrient limitations. We extend the LPJmL model to cover the terrestrial nitrogen cycle by explicitly adding processes of soil nitrogen dynamics, plant uptake, nitrogen allocation, response of photosynthesis, transpiration, and maintenance respiration to variable nitrogen concentrations in plant organs, and agricultural nitrogen management. Our implementation is based on previous model implementations (Parton et al., 2001; Gerber et al., 2010; Smith et al., 2014), but some soil processes, e.g., for denitrification and volatilization, are more complex than in Smith et al. (2014), while plant N cycling is similarly parameterized. All implemented processes and the corresponding references are described in full detail in the following sections. LPJmL is the only dynamic global vegetation model that explicitly covers natural vegetation, managed croplands and grasslands, and the full terrestrial hydrology in one consistent modeling framework (Schaphoff et al., 2018b, a). We describe all new model features in full detail and present the results of a global simulation of the historic past (1901–2009) that we use to evaluate model performance.

2 Model description

The model description focuses on the nitrogen-dependent (N-dependent) part of the model. A general description of the LPJmL model is supplied by Sitch et al. (2003), Bondeau et al. (2007), and Schaphoff et al. (2018b, a). Note that Schaphoff et al. (2018b, a) provide the most comprehensive model description available, which includes a few model features added to the model after the development of the N modules had begun and which are thus not part of the LPJmL5.0 version described here. These include several minor amendments of the code, the updated grass allocation scheme (Rolinski et al., 2018), and the updated phenology scheme for natural vegetation (Forkel et al., 2014).

In the LPJmL model vegetation is represented by different plant functional types (PFTs) that can establish concurrently within a cell. These established PFTs share the same soil stand and compete for light, water, and nitrogen resources, while crop functional types (CFTs) are established exclusively at sowing on their own soil stand.

In the predecessor version LPJmL3.5, all organic matter pools (vegetation, soil) were represented as carbon pools. We now also implemented a corresponding N pool for each of these carbon pools and pools for inorganic reactive N forms (${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$) in the soil (Fig. 1). Nitrogen dynamics have been incorporated in other dynamical vegetation models, e.g., in LPJ-GUESS (Smith et al., 2014). In addition to LPJ-GUESS our model considers not only natural vegetation but also takes into account managed crops. Furthermore, nitrogen transformation in soils is simulated in a more sophisticated way incorporating the immobilization of nitrogen. In the following sections we describe the implementation of the plant N demand, uptake, allocation, the effects of N limitation, photosynthesis, maintenance respiration, and N inputs, and transformations and losses in and/or from soils. All processes are computed at a daily time step, except for fire events (annual) and the allocation of carbon and N in plants, which is computed daily only for crops but annually for natural vegetation and before each harvest event for managed grasslands. Soil processes are vertically resolved in six soil layers including one bedrock layer.

Figure 1 Carbon and nitrogen pools and associated processes for the example of crops.

2.1 Nitrogen demand

Daily photosynthesis and maximum carboxylation capacity (V_max) are computed based on absorbed photosynthetically active radiation (APAR) and canopy conductance reflecting the level of water stress (Sitch et al., 2003). This water-stressed carboxylation capacity V_max determines the demand for N of trees, grasses, and crops in the leaves. Depending on PFT-specific requirements for V_max, the N demand of leaf, N_leaf (g N m⁻²), is calculated according to Smith et al. (2014) based on Haxeltine and Prentice (1996) as

$$\begin{array}{ll}\text{(1)}& {\displaystyle }{\text{N}}_{\mathrm{leaf}}& {\displaystyle }=\mathrm{25}\cdot \mathrm{0.02314815}/\mathrm{daylength}\cdot V_{max}{\displaystyle }& {\displaystyle }\cdot \exp (-\mathrm{0.02}\cdot (T-\mathrm{25}\left)\right)\cdot f_{\mathrm{LAI}}\left(\mathrm{LAI}\right)+\mathrm{0.00715}\cdot {\text{C}}_{\mathrm{leaf}},\end{array}$$

where C_leaf is the actual leaf carbon content (g C m⁻²) and daylength is the duration of daylight (h). The function f_LAI(LAI) is a modifier dependent on current leaf area index (LAI) accounting for a stronger leaf N content decline with canopy depth compared to incoming sunlight.

$$\begin{array}{}\text{(2)}& f_{\mathrm{LAI}}\left(\mathrm{LAI}\right)=\left\{\begin{array}{ll}max(\mathrm{0.1},\mathrm{LAI})& \text{for}\mathrm{LAI}<\mathrm{1}\\ \exp (\mathrm{0.08}\cdot min(\mathrm{LAI},\mathrm{7}\left)\right)& \text{otherwise}\end{array}\right.\end{array}$$

The pre-factor 0.12 in the exponential term of Smith et al. (2014) has been replaced by 0.08 for two reasons. First, we find that canopy C : N ratios are too low for the original value. Second, the computed values for the average leaf C : N ratio of the canopy should monotonically increase with LAI, whereas they decline again at higher LAI. This unwanted decline is not completely prevented with our pre-factor of 0.08 but much weaker and occurs only at much higher LAI values than in the original implementation (see Fig. S1 in the Supplement). We choose a maximum of LAI=7 and for LAI<1 a linear decrease to avoid too-high respiration rates at low LAI levels, where C : N ratios would become very small otherwise. Daily gross photosynthesis A_gd depends on the light-limited photosynthesis rate J_E and Rubisco-limited photosynthesis rate J_C:

$$\begin{array}{}\text{(3)}& A_{\mathrm{gd}}={\displaystyle \frac{\left(J_E+J_C-\sqrt{(J_E+J_C{)}^{\mathrm{2}}-\mathrm{4}\cdot \mathit{\theta }\cdot J_E\cdot J_C}\right)}{(\mathrm{2}\cdot \mathit{\theta })\cdot \mathrm{daylength}}},\end{array}$$

where θ is the shape parameter describing the co-limitation of light and Rubisco activity. The value of θ of LPJmL3.5 has been changed from 0.7 to 0.9, which is in better agreement with Collatz et al. (1990) and results in lower Rubisco demand to reach the light-limited photosynthesis rate (see Fig. S2). The factor α_a determining the fraction of photosynthetic active radiation (PAR) assimilated at ecosystem level relative to leaf level has been changed from 0.5 to 0.6 to counterbalance the reduction of GPP due to additional nitrogen limitations.

Because the allocation of carbon and nitrogen for grass and tree PFTs is done on a yearly time interval, the actual carbon stored in leaves C_leaf,t at time t of the current year is calculated from the carbon stored in leaves at the end of the previous year C_leaf:

$$\begin{array}{}\text{(4)}& {\text{C}}_{\mathrm{leaf},t}={\text{C}}_{\mathrm{leaf}}+f_{\mathrm{leaf}}\cdot \sum _{t^{\prime }=\mathrm{1}}^t{\mathrm{NPP}}_{t^{\prime }},\end{array}$$

where ${\sum }_{t^{\prime }}{\mathrm{NPP}}_{t^{\prime }}$ is the accumulated biomass increment and f_leaf is the fraction of biomass that was allocated to leaves at the end of the previous year. Then the total N demand is determined by the actual (t) carboxylation-based demand for N in leaves N_leaf,t (see Eq. 1), the current N content of the other organs (roots N_root and sapwood N_sapwood for trees), and the approximated N demand for the newly accumulated NPP (Eq. 5). For this approximation, we use the allocation shares of the previous year (f_root, f_sapwood):

$$\begin{array}{ll}\text{(5)}& {\displaystyle }{\text{N}}_{\mathrm{demand},t}=& {\displaystyle }{\text{N}}_{\mathrm{leaf},t}+{\text{N}}_{\mathrm{root}}+{\text{N}}_{\mathrm{sapwood}}+{\displaystyle \frac{{\text{N}}_{\mathrm{leaf}}}{{\text{C}}_{\mathrm{leaf}}}}{\displaystyle }& {\displaystyle }\cdot (f_{\mathrm{root}}/R_{\mathrm{1}}+f_{\mathrm{sapwood}}/R_{\mathrm{2}})\cdot \sum _{t^{\prime }=\mathrm{1}}^t{\mathrm{NPP}}_{t^{\prime }},\end{array}$$

where R₁ and R₂ are the prescribed PFT-specific C : N ratios of roots and sapwood relative to the C : N ratio of leaves (Table 1).

Table 1 C : N ratios relative to the leaf C : N ratio R_i for the different plant compartments.

The daily allocation scheme of crops enables the calculation of nitrogen demand by using the carbon compartment itself. Plants maintain a store of labile N, N_store (g N m⁻²), to buffer fluctuations between N demand and supply from the soil mineral N pool (Smith et al., 2014). N demand is therefore increased by a factor of k_store=1.15 for trees and k_store=1.3 for grass and crops. Thus, the optimum N uptake fulfilling the demand N_uptake,opt can be calculated from the demand increment.

$$\begin{array}{}\text{(6)}& {\text{N}}_{\mathrm{uptake},\mathrm{opt}}=({\text{N}}_{\mathrm{demand},t}-{\text{N}}_{\mathrm{demand},t-\mathrm{1}})\cdot k_{\mathrm{store}}\end{array}$$

2.2 Nitrogen uptake

The mechanism for the uptake of N (N_uptake in g N m⁻² d⁻¹) is the same for trees, crops, and grasses. Following Smith et al. (2014), plant N_uptake is determined by soil mineral N concentrations, fine root mass, soil temperature and porosity, and plant demand for N. This is computed for all soil layers individually and summed up to compute overall N uptake:

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

where N_up,root is the maximum N uptake rate per unit of fine root mass in each layer, f_N(N_avail) parameterizes the dependence on available N, f_T(T_soil) parameterizes the temperature dependence, f_NC parameterizes the dependence on plant N : C ratio, C_root is the carbon stored in the roots, n_soillayer is the number of soil layers (n_soillayer=6), and rootdist_l determines the fraction of roots in each layer. N_up,root is $\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{3}}$ g N g C⁻¹ d⁻¹ for trees and $\mathrm{5.51}\times {\mathrm{10}}^{-\mathrm{3}}$ g N g C⁻¹ d⁻¹ for crops and grasses (Smith et al., 2014). The available N is the sum of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the soil layer l.

$$\begin{array}{}\text{(8)}& {\text{N}}_{\mathrm{avail},l}={\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}^{-}+{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l}^{+}\end{array}$$

The function f_N can be parameterized by Michaelis–Menten kinetics:

$$\begin{array}{}\text{(9)}& f_{\text{N}}\left({\text{N}}_{\mathrm{avail},l}\right)=k_{\text{N},min}+{\displaystyle \frac{{\text{N}}_{\mathrm{avail},\mathrm{l}}}{{\text{N}}_{\mathrm{avail},l}+K_{\text{N},min}\cdot {\mathit{\theta }}_{max}\cdot d_{\mathrm{soil},l}}},\end{array}$$

where d_soil,l is the soil column depth (m), θ_max is the soil-type-specific fractional pore space (dimensionless), $K_{\text{N},min}$ is 1.48 g N m⁻³ for woody and 1.19 for grassy PFTs (half-saturation concentration of fine root N uptake), and $k_{\text{N},min}$ (dimensionless) is 0.05, which is the basal rate of N uptake that is not associated with Michaelis–Menten kinetics. The function f_NC(NC_plant) is from Zaehle et al. (2010b):

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

where NC_leaf,low and NC_leaf,high are the lower and upper limits of N : C ratios and NC_plant is the actual plant N : C ratio. The lower and upper limits NC_leaf,low and NC_leaf,high are derived from the TRY database (Kattge et al., 2011). Their reciprocal C : N values for each PFT are shown in Table 2. The actual plant N : C ratio is calculated according to

$$\begin{array}{}\text{(11)}& {\mathrm{NC}}_{\mathrm{plant}}={\displaystyle \frac{{\text{N}}_{\mathrm{leaf}}+{\text{N}}_{\mathrm{root}}}{{\text{C}}_{\mathrm{leaf}}+{\text{C}}_{\mathrm{root}}}}.\end{array}$$

The temperature function f_T for N uptake is given by Thornley (1991):

$$\begin{array}{}\text{(12)}& f_T\left(T_{\mathrm{soil},l}\right)={\displaystyle \frac{(T_{\mathrm{soil},l}-T_{\mathrm{0}})\cdot (\mathrm{2}\cdot T_m-T_{\mathrm{0}}-T_{\mathrm{soil},l})}{(T_r-T_{\mathrm{0}})\cdot (\mathrm{2}\cdot T_m-T_{\mathrm{0}}-T_r)}},\end{array}$$

where $T_{\mathrm{0}}<T_r<\mathrm{2}\cdot T_m-T_{\mathrm{0}}$. For the chosen T_m=15 ^∘C, T_r=15 ^∘C, and $T_{\mathrm{0}}=-\mathrm{25}$ ^∘C, the maximum of 1 is reached at 15^∘ and the function is positive above −25 ^∘C.

The root distribution rootdist_l can be calculated from the proportion of roots from the surface to soil depth z, rootdist_z, as in Jackson RB et al. (1996):

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

where β_root is a PFT-specific parameter (for parameter values, see Table 2); rootdist_l is then given by the difference ${\mathrm{rootdist}}_{z\left(l\right)}-{\mathrm{rootdist}}_{z(l-\mathrm{1})}$. If the soil depth of the layer l is greater than the thawing depth then rootdist_l is reduced accordingly. The nonzero rootdist_l values are rescaled so that their sum is normalized to 1, accounting for the modified root distribution under freezing conditions. Soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and soil ${\mathrm{NH}}_{\mathrm{3}}^{-}$ pools are reduced accordingly every simulation day t.

$$\begin{array}{ll}\text{(14)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t+\mathrm{1}}^{-}={\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}\cdot \left(\mathrm{1}-{\mathrm{rootdist}}_l\cdot {\displaystyle \frac{N_{\mathrm{uptake}}}{{\sum }_{l=\mathrm{1}}^{n_{\mathrm{soillayer}}}{N}_{\mathrm{avail},l}}}\right)\\ \text{(15)}& {\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t+\mathrm{1}}^{+}={\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4},{\mathrm{soil}}_l,t}^{+}\cdot \left(\mathrm{1}-{\mathrm{rootdist}}_l\cdot {\displaystyle \frac{N_{\mathrm{uptake}}}{{\sum }_{l=\mathrm{1}}^{n_{\mathrm{soillayer}}}{N}_{\mathrm{avail},l}}}\right)\end{array}$$

Table 2 PFT-specific β_root based on Schaphoff et al. (2018b) and minimum and maximum leaf C : N ratios based on the TRY database (Kattge et al., 2011) with data from Kurokawa and Nakashizuka (2008), Garnier et al. (2007), Penuelas et al. (2010a, b), Fyllas et al. (2009), Loveys et al. (2003), Han et al. (2005), Ordonez et al. (2010), Atkin et al. (1999), White et al. (2000), Xu and Baldocchi (2003), Freschet et al. (2010a, b), Laughlin et al. (2010), Niinemets (2001, 1999), Willis et al. (2010), Baker et al. (2009), Patiño et al. (2009), Pakeman et al. (2009, 2008), Fortunel et al. (2009), Cornelissen et al. (1996, 1997, 2004), Quested et al. (2003), Sardans et al. (2008b, a), and Ogaya and Penuelas (2003, 2006, 2007, 2008). The C : N ratios for C₃ and C₄ grasses and crops are based on White et al. (2000).

2.3 Determination of the N limitation scalar

For trees, grass, and crops, the N limitation scalar v_scal is calculated as the ratio of N demand N_uptake,opt to actual N uptake:

$$\begin{array}{}\text{(16)}& v_{\mathrm{scal}}=min({\text{N}}_{\mathrm{uptake}}/{\text{N}}_{\mathrm{uptake},\mathrm{opt}},\mathrm{1}).\end{array}$$

The scalar v_scal is used to account for N limitation in the allocation of N to different plant organs (Sect. 2.5) and is computed as the growing season mean, which is re-initialized to zero every year for natural vegetation and at sowing for crops.

2.4 Photosynthesis and gross and net primary production under N limitation

To calculate the limitation by N availability, N stress is calculated after determining water stress on photosynthesis. If N demand from the water-limited photosynthesis rate cannot be fulfilled by N uptake, carboxylation capacity V_max has to be reduced. The reduced V_max is determined by solving Eq. (1) for V_max. Water demand is then recalculated using the reduced V_max. From this reduced V_max, the actual photosynthesis rate and canopy conductance can be calculated (Fig. 2). For the determination of the canopy conductance we assume higher PFT-specific minimum canopy conductances g_min (see Table S1 in the Supplement) than Schaphoff et al. (2018b), which are in the range of values reported by Barnard and Bauerle (2013). Furthermore, we have adjusted some additional parameters (Tables S1, S2) to meet global and local evapotranspiration fluxes under nitrogen limitation effects on transpiration.

Figure 2 Calculation of N stress of plants.

The gross primary production (GPP) derived from the actual photosynthesis rate is reduced by leaf, root, and sapwood (for tree PFTs) respiration R_leaf, R_root, and R_sapwood in order to get the net primary production (NPP). Respiration rates of roots and sapwood are assumed to be linearly dependent on the N : C ratio of the corresponding pool, whereas the respiration rate of leaves (R_leaf) is a fraction (1.5 % for C₃ plants, 3.5 % for C₄ plants) of V_max (Sitch et al., 2003):

$$\begin{array}{}\text{(17)}& {\displaystyle }{R}_{\mathrm{root}}& {\displaystyle }=k_{\mathrm{resp}}\left(T_{\mathrm{soil}}\right)\cdot {\text{N}}_{\mathrm{root}},\text{(18)}& {\displaystyle }{R}_{\mathrm{sapwood}}& {\displaystyle }=k_{\mathrm{resp}}\left(T_{\mathrm{air}}\right)\cdot {\text{N}}_{\mathrm{sapwood}},\end{array}$$

where k_resp(T) is a temperature-dependent respiration rate (g C g N⁻¹ d⁻¹) (as in Sitch et al., 2003). Therefore, higher N : C ratios lead to a reduction in net primary production (NPP), which is computed as

$$\begin{array}{}\text{(19)}& \mathrm{NPP}=\mathrm{GPP}-R_{\mathrm{growth}}-R_{\mathrm{leaf}}-R_{\mathrm{root}}-R_{\mathrm{sapwood}},\end{array}$$

where R_growth is 25 % of GPP and R_sapwood is zero for all non-woody plants.

2.5 Nitrogen allocation and turnover in plants

Carbon allocation to plant compartments follows functional and allometric rules as described by Sitch et al. (2003) and is computed annually for natural vegetation and daily for crops (Bondeau et al., 2007). The allocation rules account for the functional relationships that leaf area needs to be supported by sufficient sapwood (in trees) and fine root biomass. Fine root biomass increases relative to leaf biomass under water stress and also under nitrogen limitation. The allometric rules specify the relationship of stem diameter to plant height and crown diameter (Sitch et al., 2003). Plants require N in varying amounts to satisfy organ-specific C : N ratios. Leaf N content is determined by photosynthetic potential and structural requirements and can vary within PFT-specific limits of C : N ratios. The PFT-specific range of possible C : N ratios is based on the TRY database (Kattge et al., 2011, Table 2).

The allocation of N (N_inc) to plant compartments follows the allocation rules for carbon and ensures distribution between plant compartments as established with the relative ratios given for the C : N ratio of, e.g., roots in comparison to leaves (CN_root / CN_leaf). These relative ratios for natural vegetation are taken from Friend et al. (1997, Table 4).

For crops the C : N ratios for the storage organ are derived from Bodirsky et al. (2012). Therefore, average crop-functional-type-specific leaf C : N ratios as simulated by LPJmL5.0 were used to estimate the factors R₃ that relate leaf C : N ratios to storage organ C : N ratios (Table 1).

The allocation scheme follows the algebraic solution of the following set of equations when there are n plant compartments.

$$\begin{array}{ll}\text{(20a)}& {\displaystyle \frac{{\text{N}}_{\mathrm{1}}+a_{\mathrm{1}}\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_{\mathrm{1}}}}& {\displaystyle }=R_{\mathrm{1}}\cdot {\displaystyle \frac{{\text{N}}_{\mathrm{2}}+a_{\mathrm{2}}\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_{\mathrm{2}}}}\text{(20b)}& {\displaystyle \frac{{\text{N}}_{\mathrm{1}}+a_{\mathrm{1}}\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_{\mathrm{1}}}}& {\displaystyle }=R_{\mathrm{2}}\cdot {\displaystyle \frac{{\text{N}}_{\mathrm{3}}+a_{\mathrm{3}}\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_{\mathrm{3}}}}{\displaystyle }& {\displaystyle }\mathrm{⋮}\\ \text{(20c)}& {\displaystyle \frac{{\text{N}}_{\mathrm{1}}+a_{\mathrm{1}}\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_{\mathrm{1}}}}& {\displaystyle }=R_{n-\mathrm{1}}\cdot {\displaystyle \frac{{\text{N}}_n+a_n\cdot {\text{N}}_{\mathrm{inc}}}{{\text{C}}_n}}\text{(20d)}& {\displaystyle }\sum _{i=\mathrm{1}}^n{a}_i& {\displaystyle }=\mathrm{1}\end{array}$$

${\text{C}}_{\mathrm{1}},{\text{C}}_{\mathrm{2}},\mathrm{\dots },{\text{C}}_n$, N₁, N₂, …, N_n are the C and N pools of plant compartments $\mathrm{1},\mathrm{2},\mathrm{\dots },n$, and $R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots }{R}_{n-\mathrm{1}}$ are the relative C : N ratios in comparison to leaves. The system is solved for $a_{\mathrm{1}},a_{\mathrm{2}},\mathrm{\dots }{a}_n$ so that the relative ratios $R_{\mathrm{1}},\mathrm{\dots },R_{n-\mathrm{1}}$ are ensured. Thus, the model has to solve the equation system for n=2 pools for grass, n=3 pools for trees, and n=4 pools for crops. If the N : C ratio for a pool is below the PFT-specific minimum N : C ratio allowed, then the excess carbon is put into the litter pools. To avoid overly large C fluxes from excess carbon to the litter pools in N-limited environments, we have introduced a sink limitation for the photosynthesis of trees. For this, the excess carbon from the sapwood pool is stored in an additional carbon pool C_excess. If this excess pool is filled and if there is a minimum C_sapwood pool of at least 1 kg m⁻², photosynthesis is downregulated by a scaling factor s in the following year (Eq. 21). At the end of the year, the newly acquired carbon (NPP) and the C_excess are allocated to the plant organs according to the usual allocation rules. If all carbon can be allocated within allowed compartment-specific C : N ratios, the C_excess pool is empty afterwards and photosynthesis is no longer downregulated.

$$\begin{array}{ll}\text{(21)}& {\displaystyle }s& {\displaystyle }=(\mathrm{1}+K_{\text{M}})\cdot {\displaystyle \frac{f}{f+K_{\text{M}}}},f{\displaystyle }& {\displaystyle }=min\left(\mathrm{1},{\displaystyle \frac{N_{\mathrm{sapwood}}}{{\text{C}}_{\mathrm{sapwood}}+{\text{C}}_{\mathrm{excess}}}}\cdot {\displaystyle \frac{R_{\mathrm{2}}}{{\mathrm{NC}}_{\mathrm{leaf},\mathrm{low}}}}\right),\end{array}$$

where K_M=0.1 is the Michaelis constant of the Michaelis–Menten kinetics and R₂ is the relative C : N ratio of sapwood with respect to leaves.

Similar to water stress, we assume that plants allocate more biomass to roots under N limitation. For this, the leaf to root mass ratio (lmtorm) is modified by the minimum of the N limitation factor v_scal and the water limitation factor w_scal. Both factors are computed as growing season means with daily updates, i.e., for the entire calendar year for natural vegetation, between harvest events for managed grasslands, and since sowing for crops.

LPJmL employs PFT-specific turnover rates for living leaves and fine roots. At turnover the corresponding amount of carbon is moved into the litter pools, whereas not all of the associated N is disposed of but remains in the plant. We assume that grasses and deciduous trees recover k_turn=70 % of their N upon biomass turnover, whereas evergreen trees only recover k_turn=20 %. At turnover sapwood carbon is transformed into heartwood carbon. Not all nitrogen from sapwood turnover goes into heartwood, and only a fraction f_heartwood=0.7 of nitrogen is transformed.

2.6 Nitrogen transformation in soils

Nitrogen occurs in soils in different reactive forms, mainly the organic forms nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), which are represented by different pools in LPJmL5.0. Transformations between different forms of N in the soil are represented by mineralization, immobilization, nitrification, and denitrification and are simulated in sequential order. Each soil and litter pool consists of carbon and nitrogen stocks and the resulting C : N ratios are flexible. Losses from the soil are represented by the implemented nitrification, leaching, denitrification, and volatilization processes. The corresponding pools and fluxes are depicted in Fig. 3 and described, including their parameterization (see Table S2), in this section.

Figure 3 Nitrogen transformations and losses in soils. Pools and fluxes are denoted by boxes and arrows, respectively.

2.6.1 Mineralization of nitrogen

The mineralization of N from soil organic matter and the decomposition of litter pools follow that of carbon as described by Schaphoff et al. (2013). First, for each soil layer the fluxes of carbon from the soil into the atmosphere are calculated and the respective fluxes of N, reflecting the actual C : N ratios of the material, are transferred to the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ soil pool of the corresponding soil layer.

Fluxes (F) of carbon and nitrogen for slow (s) and fast (f) pools (P) depend on parameters $k_{\mathrm{soil}\mathrm{10}}^f=\mathrm{0.03}$, $k_{\mathrm{soil}\mathrm{10}}^s=\mathrm{0.001}$ (per year), and R(T,M) as a function of temperature (T) and soil moisture (M) per soil layer (l).

$$\begin{array}{ll}\text{(22)}& {\displaystyle }& {\displaystyle }{F}_l^x=max(\mathrm{0},P_l^x\cdot (\mathrm{1}-\exp (-k_{\mathrm{soil}\mathrm{10}}^x\cdot R(T_l,M_l\left)\right)),{\displaystyle }& {\displaystyle }x\in (s,f),\end{array}$$

where

$$\begin{array}{ll}\text{(23)}& {\displaystyle }R(T_l,M_l)=& {\displaystyle }{T}_l\cdot (\mathrm{0.04021601}-\mathrm{5.00505434}\cdot M_l^{\mathrm{3}}{\displaystyle }& {\displaystyle }+\mathrm{4.26937932}\cdot M_l^{\mathrm{2}}+\mathrm{0.71890122}\cdot M_l).\end{array}$$

The mineralization of soil N, N${}_{\mathrm{miner},\mathrm{soil},l}$, in soil layer l is given by

$$\begin{array}{}\text{(24)}& {\text{N}}_{\mathrm{miner},\mathrm{soil},l}=F_l^f+F_l^s.\end{array}$$

Whereas the mineralization fluxes of carbon go completely to the atmosphere as CO₂, mineralized N goes to the mineral pools, where it is subject to further transformation (Parton et al., 2001).

The decomposition of N in soil organic material (N_decom) consists of a mineralization part (A_f=0.6, dimensionless) that forms ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and a humification part (1−A_f), in which organic N from the litter pool is transferred to the soil pools. The humification flux is divided into fluxes to slow (s) and fast (f) N soil pools (P) for which the parameter F_f=0.98 (dimensionless) specifies the portion that goes to the fast soil pool.

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }{P}_{l,t+\mathrm{1}}^f=P_{l,t}^f+F_f\cdot (\mathrm{1}-A_f)\cdot {\text{N}}_{\mathrm{decom}}\cdot {\text{N}}_{\mathrm{shift},l}^f,\text{(26)}& {\displaystyle }& {\displaystyle }{P}_{l,t+\mathrm{1}}^s=P_{l,t}^s+(\mathrm{1}-F_f)\cdot (\mathrm{1}-A_f)\cdot {\text{N}}_{\mathrm{decom}}\cdot {\text{N}}_{\mathrm{shift},l}^s,\end{array}$$

where the annual shift rates N${}_{\mathrm{shift},l}^{s,f}$ describe the organic matter input from the different PFTs into the respective layer due to cryoturbation and bioturbation (Schaphoff et al., 2013).

Net mineralized material N${}_{\mathrm{miner},\mathrm{litter},l}$ is

$$\begin{array}{ll}\text{(27)}& {\displaystyle }& {\displaystyle }{\text{N}}_{\mathrm{miner},\mathrm{litter},l}={\displaystyle }& {\displaystyle }{A}_f\cdot {\text{N}}_{\mathrm{decom}}\cdot \left(F_f\cdot {\text{N}}_{\mathrm{shift},l}^f+(\mathrm{1}-F_f)\cdot {\text{N}}_{\mathrm{shift},l}^s\right),\end{array}$$

which adds N to an intermediate N mineralization pool

$$\begin{array}{}\text{(28)}& {\text{N}}_{\mathrm{miner},l}={\text{N}}_{\mathrm{miner},\mathrm{soil},l}+{\text{N}}_{\mathrm{miner},\mathrm{litter},l}.\end{array}$$

In contrast to Parton et al. (2001) in which 20 % of this pool is directly nitrified to ${\mathrm{NO}}_{\mathrm{3}}^{-}$, we follow Schimel and Bennett (2004) and transfer all mineralized N to the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool.

$$\begin{array}{}\text{(29)}& {\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t+\mathrm{1}}^{+}={\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t}^{+}+{\text{N}}_{\mathrm{miner},l}\end{array}$$

2.6.2 Nitrogen immobilization

Immobilization, i.e., the transformation of mineral N to organic N in soils, is determined per soil layer directly after soil and litter mineralization, following the LM3V land model described by Gerber et al. (2010). If available mineral soil N is constraining immobilization, mineral N is first immobilized into the fast soil pool and then into the slow soil pool. The immobilized N, N_immo,l, is calculated according to

$$\begin{array}{ll}\text{(30)}& {\displaystyle }{\text{N}}_{\mathrm{immo},l}=& {\displaystyle }{F}_f\cdot (\mathrm{1}-A_f)\cdot \left({\text{C}}_{\mathrm{decom}}/{\mathrm{CN}}_{\mathrm{soil}}-N_{\mathrm{decom}}\right){\displaystyle }& {\displaystyle }\cdot N_{\mathrm{shift},l}^f\cdot {\displaystyle \frac{N_{\mathrm{sum},l}/d_{\mathrm{soil},l}}{k_{\text{N}}+N_{\mathrm{sum},l}/d_{\mathrm{soil},l}}},\end{array}$$

where CN_soil is the desired soil C : N ratio of 15 (dimensionless) for all soil types, d_soil,l is the soil depth of layer l in meters, $k_{\text{N}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ (g N m⁻³) is the half-saturation concentration for immobilization in soils (Gerber et al., 2010), and N${}_{\mathrm{shift},l}^f$ is the parameter that determines the distribution of the humified organic matter in the topsoil to the different soil layers l (Schaphoff et al., 2013). The available mineral N in the soil layer l (N_sum,l in g N m⁻²) is the sum of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$.

$$\begin{array}{}\text{(31)}& {\text{N}}_{\mathrm{sum},l}={\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l}^{+}+{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}^{-}\end{array}$$

The immobilized N (N_immo,l) is added to the fast soil N pool of layer l and subtracted from the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pools.

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{P}_{\mathrm{soil},l,t+\mathrm{1}}^f=P_{\mathrm{soil},l,t}^f+min(N_{\mathrm{immo},l},N_{\mathrm{sum},l})\end{array}$$

$$\begin{array}{ll}\text{(33)}& {\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t+\mathrm{1}}^{+}={\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t}^{+}-{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t}^{+}\cdot min(N_{\mathrm{immo},l}/N_{\mathrm{sum},l},\mathrm{1})\\ \text{(34)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t+\mathrm{1}}^{-}={\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}-{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}\cdot min(N_{\mathrm{immo},l}/N_{\mathrm{sum},l},\mathrm{1})\end{array}$$

The immobilization into the slow soil N pool ($P_{\mathrm{soil},l,t+\mathrm{1}}^s$) is computed accordingly as in Eq. (30) but with (1−F_f) instead of F_f.

2.6.3 Nitrification

Nitrogen fluxes from nitrification in the soil are modeled modified after Parton et al. (2001) with the schematic representation of a series of pipes for the main flow from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ over ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to N₂ from which N₂O leaks in between. As suggested by Parton et al. (2001, Eq. 2), nitrification is computed as a fixed fraction of the mineralization flux (see Sect. 2.6.1) and an explicit transformation flux $F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ from ammonium to nitrate in g N m⁻² d⁻¹, which is described here.

$$\begin{array}{}\text{(35)}& F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}=K_{max}\cdot F_{\mathrm{1}}\left(T_{\mathrm{soil},l}\right)\cdot F_{\mathrm{1}}\left(W_{\mathrm{sat},l}\right)\cdot F\left(\mathrm{pH}\right)\cdot {\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l}^{+},\end{array}$$

where ${\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l}^{+}$ is the model-derived soil ammonium concentration (g N m⁻²), K_max is the maximum nitrification rate of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ ($K_{max}=\mathrm{0.1}$ d⁻¹), F₁(T_soil,l) is the limiting function for temperature, and F₁(W_sat,l) the corresponding limiting function for water saturation W_sat,l. Parton et al. (2001) show nitrification rates after data from Malhi and McGill (1982) in Table 3 without a formula. Using these data from three different sites in the US, Canada, and Australia, we fitted a bell-shaped function for the temperature dependence:

$$\begin{array}{}\text{(36)}& F_{\mathrm{1}}\left(T_{\mathrm{soil},l}\right)=\exp (-(T_{\mathrm{soil},l}-a{)}^{\mathrm{2}}/(\mathrm{2}\cdot b^{\mathrm{2}})),\end{array}$$

where a=18.79 ^∘C and b=5.26 give the best fist to the data (see Fig. S3). The function is also applicable for negative values.

The soil water response function F₁(W_sat) is parameterized according to Doran et al. (1988) as described in Parton et al. (1996):

$$\begin{array}{ll}\text{(37)}& {\displaystyle }{F}_{\mathrm{1}}\left(W_{\mathrm{sat},l}\right)=& {\displaystyle }{\left({\displaystyle \frac{W_{\mathrm{sat},l}-b_{\mathrm{nit}}}{a_{\mathrm{nit}}-b_{\mathrm{nit}}}}\right)}^{d_{\mathrm{nit}}\cdot (b_{\mathrm{nit}}-a_{\mathrm{nit}})/(a_{\mathrm{nit}}-c_{\mathrm{nit}})}{\displaystyle }\cdot & {\displaystyle }{\left({\displaystyle \frac{W_{\mathrm{sat},l}-c_{\mathrm{nit}}}{a_{\mathrm{nit}}-c_{\mathrm{nit}}}}\right)}^{d_{\mathrm{nit}}},\end{array}$$

where W_sat,l is the water-filled pore space of soil layer l, and parameters a_nit to d_nit are given for sandy and medium soil (Table S2).

This soil pH function is based on Parton et al. (1996).

$$\begin{array}{}\text{(38)}& F\left(\mathrm{pH}\right)=\mathrm{0.56}+\arctan (\mathit{\pi }\cdot \mathrm{0.45}\cdot (-\mathrm{5}+\mathrm{pH}\left)\right)/\mathit{\pi }\end{array}$$

Soil pH values are taken from the WISE data set (Batjes, 2000). Part of the N during nitrification is lost to the atmosphere as nitrous oxide (N₂O). Parton et al. (2001) assume that the N₂O flux $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (in g N m⁻² d⁻¹) is proportional to the nitrification rate with

$$\begin{array}{}\text{(39)}& F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=K_{\mathrm{2}}\cdot F_{{\mathrm{NO}}_{\mathrm{3}}^{-}},\end{array}$$

where K₂ is the fraction of nitrified N lost as N₂O flux (K₂=0.02). Finally, soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ are updated accordingly.

$$\begin{array}{}\text{(40)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t+\mathrm{1}}^{-}={\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}+(\mathrm{1}-K_{\mathrm{2}})\cdot F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\text{(41)}& {\displaystyle }& {\displaystyle }{\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t+\mathrm{1}}^{+}={\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l,t}^{+}-F_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\end{array}$$

2.6.4 Denitrification

The reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to NO₂ and N₂ is determined for each soil layer using the implementation in SWIM (Krysanova and Wechsung, 2000).

$$\begin{array}{}\text{(42)}& D_{{\mathrm{NO}}_{\mathrm{3}}^{-}}=F_{\mathrm{2}}\left(W_{\mathrm{sat},l}\right)\cdot F_{\mathrm{2}}(T_{\mathrm{soil},l},{\text{C}}_{\mathrm{org},l})\cdot {\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}^{-},\end{array}$$

where F₂(W_sat,l) is the water response function and F₂(T,C) the soil temperature and carbon reaction function. The water response function depends on the water-filled pore space W_sat,l in the following way.

$$\begin{array}{}\text{(43)}& F_{\mathrm{2}}\left(W_{\mathrm{sat},l}\right)=\mathrm{6.664096}\times {\mathrm{10}}^{-\mathrm{10}}\cdot \exp (\mathrm{21.12912}\cdot W_{\mathrm{sat},l})\end{array}$$

The water response function shows a qualitatively similar behavior to Eq. (151) from SWIM while ensuring continuity (see Fig. S4). Parameters are fitted and adjusted so that for full soil water saturation, the value is not greater than 1. The soil temperature and carbon reaction function is parameterized according to

$$\begin{array}{}\text{(44)}& F_{\mathrm{2}}(T_{\mathrm{soil},l},{\text{C}}_{\mathrm{org},l})=\mathrm{1}-\exp (-\mathrm{CDN}\cdot F_{\mathrm{2}}(T_{\mathrm{soil},l})\cdot {\text{C}}_{\mathrm{org},l}),\end{array}$$

where CDN=1.4 is the shape coefficient (Arnold et al., 2012), C_org,l is the sum of the fast and slow C pools, and F₂(T_soil,l) is the soil temperature reaction function. F₂(T_soil,l) is replaced by Eq. (C5) from Smith et al. (2014), which is only valid for positive T_soil,l. The original function from the soil and water assessment tool (SWAT) approaches 1 for high temperatures, whereas the function from Smith declines, which seems more sensible. Equation (C5) of Smith et al. (2014) is taken from Comins and McMurtrie (1993).

$$\begin{array}{ll}\text{(45)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{2}}\left(T_{\mathrm{soil},l}\right)={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{0.0326}& \text{for}{T}_{\mathrm{soil},l}\le \mathrm{0}{}^{\circ }\mathrm{C}\\ \mathrm{0.0326}+\mathrm{0.00351}\cdot T_{\mathrm{soil},l}^{\mathrm{1.652}}& \text{for}\mathrm{0}{}^{\circ }\mathrm{C}<T_{\mathrm{soil},l}<\mathrm{45.9}{}^{\circ }\mathrm{C}\\ -{\left({\displaystyle \frac{T_{\mathrm{soil},l}}{\mathrm{41.748}}}\right)}^{\mathrm{7.19}}& \\ \mathrm{0}& \text{for}{T}_{\mathrm{soil},l}\ge \mathrm{45.9}{}^{\circ }\mathrm{C}\end{array}\right.\end{array}$$

Bessou et al. (2010) assume that the N₂O flux from ${\mathrm{NO}}_{\mathrm{3}}^{-}$, $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (in g N m⁻² d⁻¹), is proportional to the denitrification rate $D_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ with

$$\begin{array}{}\text{(46)}& F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=r_{\mathrm{mx}}\cdot D_{{\mathrm{NO}}_{\mathrm{3}}^{-}},\end{array}$$

where r_mx=0.11 is the fraction of denitrified N lost as N₂O flux. The N₂ flux $F_{{\mathrm{N}}_{\mathrm{2}}}$ is then derived by

$$\begin{array}{}\text{(47)}& F_{{\mathrm{N}}_{\mathrm{2}}}=(\mathrm{1}-r_{\mathrm{mx}})\cdot D_{{\mathrm{NO}}_{\mathrm{3}}^{-}}.\end{array}$$

The soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pools have to be reduced by the denitrification flux.

$$\begin{array}{}\text{(48)}& {\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t+\mathrm{1}}^{-}={\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}-D_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\end{array}$$

2.6.5 Nitrogen leaching and movement

Nitrate movement with water fluxes is simulated as in SWAT (Neitsch et al., 2002, 2005). Nitrate is assumed to be fully dissolved in water and moves with surface runoff, lateral runoff, and percolation water. To compute the amount of nitrate transported with the water from a soil layer, we first calculate the concentration of nitrate in the mobile water. This concentration is then multiplied by the volume of surface runoff, lateral runoff, or percolation water between soil layers or into the aquifer. The amount of nitrate leached depends on the climatic and soil conditions and on the type and intensity of soil management (e.g., plant cover, soil treatment, fertilization).

The concentration of nitrate in the mobile water ${\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile},l}$ in layer l (kg N m⁻³) is

$$\begin{array}{}\text{(49)}& {\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile},l}={\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}^{-}\cdot \left(\mathrm{1}-\exp \left(\frac{-w_{\mathrm{mobile},l}}{(\mathrm{1}-\mathit{\theta })\cdot {\mathrm{SAT}}_l}\right)\right)}{w_{\mathrm{mobile},l}}},\end{array}$$

where ${\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l}^{-}$ is the content of nitrate in layer l (g N m⁻²), w_mobile is the amount of mobile water in the layer (mm), θ=0.4 is the fraction of porosity from which anions are excluded (0.5 in Neitsch et al., 2002), and SAT_l is the saturated water content of the soil layer (mm).

The mobile water w_mobile,l in the layer l is the amount of water lost by surface runoff, lateral flow, and percolation:

$$\begin{array}{}\text{(50)}& w_{\mathrm{mobile},l}=\left\{\begin{array}{ll}{Q}_{\mathrm{surf}}+Q_{\mathrm{lat},l=\mathrm{1}}+w_{\mathrm{perc},l=\mathrm{1}}& \text{for}l=\mathrm{1}\\ Q_{\mathrm{lat},l}+w_{\mathrm{perc},l}& \text{for}l>\mathrm{1}\end{array}\right.,\end{array}$$

where Q_surf is the surface runoff (only in the topsoil layer; mm), Q_lat,l is the water discharged from the layer by lateral flow (mm), and w_perc,l is the amount of water percolating to the underlying soil layer on a given day.

Finally, the amount of nitrate that is removed with surface runoff ${{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{surf}}$ and lateral flow ${{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{lat},l}$ is calculated as

$$\begin{array}{}\text{(51)}& {\displaystyle }& {\displaystyle }{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{surf}}={\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\cdot {\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile}}\cdot Q_{\mathrm{surf}},\text{(52)}& {\displaystyle }& {\displaystyle }{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{lat},l=\mathrm{1}}={\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\cdot {\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile}}\cdot Q_{\mathrm{lat},l=\mathrm{1}},\end{array}$$

for the top layer and

$$\begin{array}{}\text{(53)}& {{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{lat},l}={\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile},l}\cdot Q_{\mathrm{lat},l}\end{array}$$

for the lower soil layers, where ${\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ is the nitrate percolation coefficient. It controls the amount of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ removed from the surface layer in runoff relative to the amount removed via percolation (Neitsch et al., 2002). The value for ${\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ can range from 0.01 to 1.0. For ${\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}\to \mathrm{0}$, the concentration of nitrate in the runoff approaches 0. For ${\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}=\mathrm{1.0}$, surface runoff has the same concentration of nitrate as the percolating water. We choose for ${\mathit{\beta }}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$ a value of 0.4.

Nitrate moved to the lower soil layer with percolation ${{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l}$ is calculated as

$$\begin{array}{}\text{(54)}& {{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l}={\mathrm{conc}}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{mobile}}\cdot w_{\mathrm{perc},l}.\end{array}$$

${{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l}$ is subtracted from the current ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the soil layer and added to the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool of the following soil layer.

$$\begin{array}{ll}\text{(55)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t+\mathrm{1}}^{-}={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}-{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l}-{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{surf}}-{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{lat},l}& \text{for}l=\mathrm{1}\\ {\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l,t}^{-}+{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l-\mathrm{1}}-{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{perc},l}-{{\mathrm{NO}}_{\mathrm{3}}^{-}}_{\mathrm{lat},l}& \text{for}l>\mathrm{1}\end{array}\right.\end{array}$$

2.6.6 Nitrogen volatilization

The volatilization of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is parameterized according to Montes et al. (2009). A convective mass transfer model is applied in which the flux varies with air temperature, air velocity over the surface, and the NH₃ concentration gradient between the ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) in solution and in the air:

$$\begin{array}{}\text{(56)}& J_{{\mathrm{NH}}_{\mathrm{3}}}=h_{\text{m}}\cdot \left([{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{gas}}-[{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{air}}\right),\end{array}$$

where $J_{{\mathrm{NH}}_{\mathrm{3}}}$ is the NH₃ volatilization flux (g NH₃-N m⁻² s⁻¹), h_m is the convective mass transfer coefficient (m s⁻¹), [NH₃]_gas is the concentration of gaseous NH₃ in equilibrium with dissolved NH₃ in solution (g NH₃ − N m⁻³ air), and [NH₃]_air is the concentration of NH₃ in ambient air (g NH₃ − N m⁻³ air), which is usually very small and can be neglected. The convective mass transfer coefficient h_m is a function of temperature T (in K), air velocity U (in m s⁻¹), and characteristic length L (m) of the emitting surface.

$$\begin{array}{}\text{(57)}& h_{\text{m}}=\mathrm{0.000612}\cdot U^{\mathrm{0.8}}\cdot T^{\mathrm{0.382}}\cdot L^{-\mathrm{0.2}}\end{array}$$

The concentration of gaseous NH₃ in equilibrium with the dissolved NH₃ is determined using Henry's law. The Henry's law K_h constant relates the concentration of dissolved NH₃ in water to an equilibrium concentration of NH₃ in the air.

$$\begin{array}{}\text{(58)}& K_{\text{h}}={\displaystyle \frac{[{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{gas}}}{[{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{solution}}}}\end{array}$$

The Henry's law constant K_h can be parameterized as a function of air temperature T_air (K).

$$\begin{array}{}\text{(59)}& K_{\text{h}}=K_{\text{h}}\left(T_{\mathrm{air}}\right)=(\mathrm{0.2138}/T_{\mathrm{air}})\cdot {\mathrm{10}}^{\mathrm{6.123}-\mathrm{1825}/T_{\mathrm{air}}}\end{array}$$

The fraction of total ammoniacal N present as NH₃ can be estimated using equilibrium thermodynamic principles:

$$\begin{array}{}\text{(60)}& {\displaystyle }{f}_{{\mathrm{NH}}_{\mathrm{3}}}=& {\displaystyle }{\displaystyle \frac{[{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{solution}}}{[{\mathrm{NH}}_{\mathrm{3}}{]}_{\mathrm{solution}}+[{\mathrm{NH}}_{\mathrm{4}}^{+}{]}_{\mathrm{solution}}}},\text{(61)}& {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{\left[{\mathrm{H}}^{+}\right]}{K_{\text{a}}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{{\mathrm{10}}^{-\mathrm{pH}}}{K_{\text{a}}}}},\end{array}$$

where K_a is the dissociation constant, [H⁺] is the proton concentration in solution, and pH $=-\log \left(\left[{\mathrm{H}}^{+}\right]\right)$. The dissociation constant K_a is parameterized as a function of temperature T (K).

$$\begin{array}{}\text{(62)}& K_{\text{a}}=K_{\text{a}}\left(T\right)={\mathrm{10}}^{\mathrm{0.05}-\mathrm{2788}/T}\end{array}$$

Then the volatilization flux F_vol (in g N m⁻² d⁻¹) is calculated according to

$$\begin{array}{ll}\text{(63)}& {\displaystyle }{F}_{\mathrm{vol}}=& {\displaystyle }\mathrm{86}\mathrm{400}\cdot h_{\text{m}}(U,T,L)\cdot K_{\text{h}}\left(T\right){\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{{\mathrm{10}}^{-\mathrm{pH}}}{K_{\text{a}}\left(T\right)}}}\cdot {\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l=\mathrm{1}}^{+}/d_{\mathrm{soil},l=\mathrm{1}},\end{array}$$

and soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is reduced in the top layer l=1 accordingly:

$$\begin{array}{}\text{(64)}& {\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l=\mathrm{1},t+\mathrm{1}}^{+}={\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l=\mathrm{1},t}^{+}-F_{\mathrm{vol}}.\end{array}$$

2.7 Nitrogen and fire

Fire creates emissions of N₂O and NO_x and leaves nutrient-rich ashes as well as charcoal. Following Gerber et al. (2010), the flux of N due to fire is divided between atmospheric emissions and ash introduction to the nitrate pool of the upper soil layer NO${}_{\mathrm{3},\mathrm{soil},l=\mathrm{1}}^{-}$.

$$\begin{array}{}\text{(65)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{fire}}={\text{C}}_{\mathrm{fire}}\cdot N_{\mathrm{pool}}/{\text{C}}_{\mathrm{pool}}\text{(66)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{emission}}=(\mathrm{1}-q_{\mathrm{ash}})\cdot N_{\mathrm{fire}}\text{(67)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l=\mathrm{1},t+\mathrm{1}}^{-}={\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l=\mathrm{1},t}^{-}+q_{\mathrm{ash}}\cdot N_{\mathrm{fire}}\end{array}$$

Above, q_ash=0.45 is the fraction of N going into the topsoil layer ${\mathrm{NO}}_{\mathrm{3}}^{-}$.

2.8 Biological N fixation

The biological fixation of N occurs at all stands with an exception for agricultural stands. There, it is applied only for the nodulating leguminous crops pulses and soybean. For these two crops, biological N fixation (BNF) is simply the difference between N demand and N uptake, basically first using the easily plant-available N from the soils and then fixing extra N at no extra cost. For natural vegetation and grasslands, the function from Cleveland et al. (1999) is applied depending on the 20-year average annual evapotranspiration (etp; in mm yr⁻¹). BNF (in g N m⁻² d⁻¹) is assumed to only occur if there is a minimum root biomass of 20 g C m⁻². All N fixed by BNF is assumed to enter the system as ammonium in the upper soil layer (l=1).

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

The function gives linearly increasing values that are positive for etp>7.35 and are set to zero otherwise. Note that in Zaehle et al. (2010a) this function is also cited in the Supplement but with a positive intercept that is not the original equation from Cleveland et al. (1999).

2.9 Nitrogen fertilization of crops

Fertilizer is applied at sowing and when the amount of fertilizer is larger than 5 g N m⁻², only half of the fertilizer is applied at sowing. The second application occurs when the phenological stage of the crop development (fphu) exceeds 0.4 to avoid large loss fluxes (leaching, volatilization, nitrification, denitrification) when fertilizing large amounts of N at the beginning of the season.

Nitrogen fertilizer is assumed to be ammonium nitrate (NH₄NO₃), so half of the applied rate is put into the topsoil layer nitrate pool (${\mathrm{NO}}_{\mathrm{3},\mathrm{soil},l=\mathrm{1}}^{-}$) and the other half into the topsoil layer ammonium pool (${\mathrm{NH}}_{\mathrm{4},\mathrm{soil},l=\mathrm{1}}^{+}$).

3 Model setup

For the assessment of model performance, we focus on the historic period 1901–2009. The spatial longitudinal–latitudinal resolution is $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$. We conduct six different sets of simulations, two simulations with the carbon-only predecessor model version LPJmL3.5 and four with the newly implemented nitrogen version LPJmL5.0. Both model versions are used for a standard historic simulation with dynamic land use change, referred to as LPJmL35 and LPJmL5, respectively, and for a simulation without human land use in which potential natural vegetation (PNV) is simulated on the entire land surface. These runs are referred to as LPJmL35-PNV and LPJmL5-PNV. For analyzing the current N limitation, we also conduct a simulation with dynamic land use but with unlimited N supply (LPJmL5-nL) and one with no fertilizer application (LPJmL5-nF). Unlimited N supply has been modeled by a deposition rate of 1 Kg N m⁻² d⁻¹ on every LPJ cell.

3.1 Model input

Model simulations are driven with observational monthly input data on monthly precipitation from the Global Precipitation Climatology Centre (GPCC Full Data Reanalysis version 7.0; Becker et al., 2013) and daily mean temperatures from the Climatic Research Unit (CRU TS version 3.23; University of East Anglia Climatic Research Unit and Harris, 2015; Harris et al., 2014). Radiation data, shortwave downward and net downward longwave, are provided by reanalysis data from ERA-Interim (Dee et al., 2011). Monthly precipitation is allocated to individual days of the corresponding month by deriving the number of wet days per month synthetically as suggested by New et al. (2000).

Land use input is derived from MIRCA2000 (Portmann et al., 2010) using the maximum monthly growing areas per crop and grid cell combined with the extent of areas equipped for irrigation (Siebert et al., 2015). HYDE3 (Klein Goldewijk and van Drecht, 2006) gives the relative changes in cropland and pasture extent backward to 1700. Further information is given by Fader et al. (2010).

The global data set “Simulated Topological Network” (STN-30) drainage direction map (Vorosmarty and Fekete, 2011) gives the transport directions of the river-routing scheme. We use the GRanD database (Lehner et al., 2011), which provides detailed information on water reservoirs that includes information on storage capacity, total area, and main purpose. Furthermore, information on natural lakes is obtained from Lehner and Döll (2004).

Nitrogen deposition is based on the ACCMIP database (Lamarque et al., 2013) for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ separately, which is applied daily to the corresponding mineral N pools of the topsoil layer. Dry and wet deposition are not distinguished. Soil pH data are taken from the WISE data set (Batjes, 2000). Fertilizer data are crop specific, but static in time. We use the data supplied by the Global Gridded Crop Model Intercomparison (GGCMI phase 1; Elliott et al., 2015), which is based on gridded mineral fertilizer data (Mueller et al., 2012) and manure data (Potter et al., 2010) from which 60 % are assumed to be plant available and thus included, whereas the remainder are ignored and not included in Elliott et al. (2015).

3.2 Model initialization, spin-up, and equilibration of soil

All carbon and water pools are initialized to zero except soil water, soil carbon, and soil temperatures, which are computed from a 30-year averaged climate. Then a spin-up simulation of 5000 years is performed to bring permafrost extent, vegetation patterns, and carbon stocks into dynamic equilibrium. The long spin-up time is necessary for reaching these equilibrium states in the permafrost regions (Schaphoff et al., 2013).

Soil N pools (organic and mineral) are initialized with assumptions to allow for initial vegetation growth. Organic N pools (slow, fast) and mineral N pools (${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$) are set to 10⁴ g N m⁻². After 1320 simulation years, vegetation composition is assumed to have reached an equilibrium (Schaphoff et al., 2013) and litterfall is tracked for another 3680 years to allow for estimating soil carbon and soil N stocks. Based on these estimates for carbon and N fluxes under equilibrium conditions, nitrogen and carbon pools are re-initialized following Sitch et al. (2003). Hence, all N from the initialization is removed that is not supporting plant growth because of other constraints such as water shortage (e.g., in deserts).

A second spin-up phase of 390 years is conducted for all versions, including land use change (except in the PNV runs LPJmL35-PNV and LPJmL5-PNV) by using the land use input of Fader et al. (2010) to capture the influence of historic land use change on the carbon and nitrogen pools in soil and vegetation.

Table 3 Global carbon pools (soil and vegetation carbon) and fluxes (net biome productivity NBP, net primary productivity NPP, and gross primary productivity GPP) for the six different experiments (averages over the period 2000 to 2009). The suffix -PNV denotes experiments with potential natural vegetation, -nL with unlimited N supply, and -nF without fertilizer input.

Figure 4 Relative global changes in GPP (a), NPP (b), vegetation carbon (c), and soil carbon (d). Relative changes are calculated by dividing the values by their 1901–1910 average to make the different model versions and settings comparable. The blue lines denote values for LPJmL3.5 with land use (LPJmL3.5), the light blue lines for LPJmL3.5 with natural vegetation only (LPJmL3.5-PNV), the red lines for LPJmL5.0 with land use (LPJML5), and the orange lines for LPJmL5.0 with natural vegetation only (LPJmL5-PNV).

4 Results

4.1 Carbon pools and fluxes

Simulations with LPJmL5.0 result in carbon pools, NPP, and GPP fluxes comparable to the carbon-only LPJmL3.5 version (Table 3) and show a similar temporal dynamic (Fig. 4). Net biome productivity (NBP) shows in both versions a carbon source driven by productivity on managed grassland. The residual sink is at the lower end of global estimations from Le Quéré et al. (2016), but land use and land use change emissions are too high for the LPJmL5.0 simulation. The actual vegetation carbon pool is strongly limited by current N levels and increases substantially across all ecosystems, when N limitations are lifted (LPJmL5-nL). Under actual N limitations and static current fertilizer levels (Elliott et al., 2015), global GPP is relatively stable throughout the simulation period (1901–2009; red line in Fig. 4a) as the expansion of cropland into increasingly low-input areas compensates for the increase in GPP in natural vegetation (orange line in Fig. 4a). NPP increases in the standard simulation with dynamic land use (LPJmL5), but not as strongly as for natural vegetation (compare red and orange lines in Fig. 4b). The difference in global annual NPP between simulations with natural vegetation only and dynamic land use increases significantly from 3 % (LPJmL35-PNV-LPJmL35) to 19 % (LPJmL5-PNV-LPJmL5). This indicates that agricultural land is increasingly N limited so that the C : N-ratio-dependent maintenance respiration declines and NPP increases, whereas GPP does not. This is in part because simulations are conducted with static fertilizer inputs and because land use change predominantly causes cropland expansion in countries with low fertilizer use so that the global average fertilizer use declines, causing higher N limitation on agricultural land. Land-use-driven declines in vegetation carbon over the 20th century are similar between the carbon-only LPJmL35 and the simulation with nitrogen LPJmL5 (Fig. 4c), but soil carbon stocks decline with nitrogen, whereas increases in the natural vegetation balance the land-use-change-induced losses in the carbon-only version (Fig. 4d).

When N limitations are lifted through unlimited N supply (LPJmL5-nL), GPP is mostly increased, except in very dry environments. Most limitation occurs in the boreal zone and in the tundra (Fig. 5a). The scatter plot (Fig. 5b) shows that the GPP increase mainly occurs in low to moderately productive areas. Decreases in GPP under unlimited N supply are possible where other factors are strongly limiting (e.g., water) and the higher N supply leads to higher maintenance respiration under lower tissue C : N ratios so that less biomass is available for leaves and thus less light can be intercepted.

Figure 5 Ratio of GPP under actual N limitations (LPJmL5) to unlimited N supply (LPJmL5-nL) (a); values less than 1 indicate higher GPP under unlimited N supply. The scatter plot (b) shows that GPP is increased through additional N supply mostly in low to moderately productive regions.

4.2 Nitrogen pools and fluxes

The comparison of global N values to literature estimates is divided between values including natural vegetation only and those considering land use dynamics (Table 4). Whereas several estimates exist for global N pools and fluxes under natural vegetation, those including land use dynamics are rather rare and given mostly for emissions from the soil (e.g., denitrification or N₂O). Lifting N limitation (LPJmL5-nL) results in similar responses for all global nitrogen pools (Table 4). Vegetation N and plant uptake increase substantially by a factor of 2.87 and 2.89, respectively, whereas soil N pools increase by a factor of 1.81. The omission of N fertilizers has little effect on N pools, which are dominated by natural ecosystems, but strongly affect nitrogen losses, especially leaching and volatilization fluxes (Table 4). A comparison to literature estimates is discussed further in Sect. 4.4.1.

4.3 Land use and nitrogen dynamics

The role of human land use for the limitation of plant growth by nitrogen availability is apparent when comparing simulations with land use (LPJmL5, red lines in Fig. 6) and natural vegetation only (LPJmL5-PNV, orange lines in Fig. 6). The nitrogen pool in the natural vegetation is stable during the 20th century (Fig. 6a) with some minor fluctuations and the global C : N ratio increases slightly by 3.5 % (Fig. 6b), whereas vegetation nitrogen with the inclusion of historical land use declines by more than 26 %. The predominant difference between the two simulations is the 22 % increase in losses of N by leaching under land use (Fig. 6c), which is caused by the additional fertilizer and irrigation water inputs under land use.

Figure 6 Simulated global time series of vegetation nitrogen (a), vegetation C : N ratios (b), and leaching (c) with land use (LPJmL5, red line) and potential natural vegetation (LPJmL5-PNV, orange line).

Figure 7 Observed C : N ratios of harvested crops versus simulated mean ratios for the crop PFTs (a) with N limitation (LPJmL5) and (b) without N limitation (LPJmL5-nL). The vertical error bars denote the 95 % percentile.

Table 4 Global nitrogen pools and fluxes for the four different experiments with LPJmL5 and literature estimates (averages over the period 2000 to 2009). The suffix -PNV denotes experiments with potential natural vegetation, -nL with unlimited N supply, and -nF without fertilizer input.

¹ Smith et al. (2014), ² Zaehle et al. (2010a), ³ Xu-Ri and Prentice (2008), ⁴ Lin et al. (2000), ⁵ Post et al. (1985), ⁶ Bouwman et al. (2013), ⁷ Sutton et al. (2013), ⁸ Galloway et al. (2013), ⁹ Galloway et al. (2004), ¹⁰ Vitousek et al. (2013).

The approximated relationships between leaf C : N ratios and storage organ C : N ratios based on Bodirsky et al. (2012) lead to consistent but variable C : N ratios in harvested crop organs, reflecting differences between crop types (Fig. 7). The leguminous crops (soybean, pulses) are not limited by N, as they can acquire the necessary N via biological N fixation. For these, the C : N ratios of harvested organs are typically underestimated. Under unlimited N supply, C : N ratios are typically reduced (Fig. 7b).

We find that agricultural land use and the associated fertilizer application greatly increases nitrogen pollution. Leaching (+93 %) and ammonia volatilization (+33 %) increase strongly, which is almost entirely driven by fertilizer inputs, not by land use change (compare LPJmL5 with LPJmL5-PNV and LPJmL5-nF in Table 4). In contrast, N₂O emissions only change slightly when agricultural land use is accounted for as increases in denitrification are partially compensated for by decreases in nitrification under reduced net mineralization (mineralization minus immobilization flux) of soil organic matter (Table 4). The effect of agricultural land use and fertilizer application is also clearly detectable in the spatial patterns of leaching. The ratio of LPJmL5-PNV to LPJmL5 (Fig. 8a) is mostly below 1, indicating higher leaching in agricultural areas. In natural vegetation under dry conditions ratios above 1 can also occur (Fig. 8b).

4.4 Model evaluation

To evaluate model performance, we focus here on carbon and nitrogen pools and fluxes at global and specific sites. Many estimates are also model based so that these cannot be used for model evaluation but only for putting our simulation results into context.

4.4.1 N pools and fluxes

Typically, simulated N pools and fluxes are within literature estimates (Table 4), although literature estimates are often broad, reflecting substantial uncertainty in these values. Values from other model simulations are generally not suitable for an evaluation of model results as they may be wrong (Kelley et al., 2013), and we only include them here for pools and fluxes for which no independent data are available. These model-based reference points include the vegetation N pool, soil mineral N pools, plant uptake rates, mineralization rates, and most global values assessed for potential natural vegetation (PNV). The vegetation N of the potential natural vegetation simulation, LPJmL5-PNV, is slightly below the other model-based estimates (Table 4), whereas other fluxes (e.g., plant uptake of N or mineralization) and pools (soil organic N) are within the range of reported values. For simulations with land use history, LPJmL5, a comparison with independent data is possible for most of the emissions from the soil. Our values for leaching and N₂ emissions are slightly below other estimates. For the complementary flux, N₂O emissions from denitrification, there is no other estimate, but total N₂O emissions from denitrification and nitrification are within the range of other estimates again (Table 4). Xu-Ri and Prentice (2008) provide the only study reporting global soil pools of mineral N forms, but for potential natural vegetation only and for the upper 1.5 m soil layer, and this estimate is based on their model application, not on independent data. Our values for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ are the inventory of 3 m soil. In comparison to Xu-Ri and Prentice (2008) we overestimate ${\mathrm{NO}}_{\mathrm{3}}^{-}$ values by a factor of 5 and underestimate ${\mathrm{NH}}_{\mathrm{4}}^{+}$ values by a factor of 2 in the soil. The accumulation of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in lower soil layers is a phenomenon also reported by Walvoord et al. (2003) and Ascott et al. (2017). In any case, this nitrogen pool is largely inaccessible to plants, as they have very little root access to these layers in our model. Also, higher ${\mathrm{NO}}_{\mathrm{3}}^{-}$ than ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations are common in soils as reported by, e.g., Kabala et al. (2017).

Figure 8 Ratio of leaching flux under potential natural vegetation (LPJmL5-PNV) to the leaching flux under actual land use patterns (LPJmL5) (a); values less than 1 indicate higher leaching under actual land use patterns. The scatter plot in (b) shows that leaching is increased strongly mostly in regions where leaching is low under potential natural vegetation.

Figure 9 Observed GPP (a), NPP (b), and vegetation carbon (c) data (Luyssaert et al., 2007) versus simulated data for natural vegetation (LPJmL5-PNV). The horizontal error bars denote the minima and maxima of observed site data belonging to the same LPJmL grid cell and the open dot is the mean. The thin black line denotes the 1 : 1 line. RMSE: root mean square error, NMSE: normalized mean square error, NME: normalized mean error, s: slope, p: significance.

Figure 10 Maize yield simulations (in tons of fresh matter (FM) ha⁻¹) for the 10 top-producing countries for the carbon-only LPJmL3.5 version (LPJmL35), the version with N limitation (LPJmL5), and with unlimited N supply (LPJmL5-nL). The residuals plotted are the detrended observed and simulated yields.

4.4.2 Carbon cycle dynamics

Carbon dynamics are mostly unchanged from the predecessor version LPJmL3.5 (Table 3 and Fig. 4). In comparison to measured site-level GPP, NPP, and vegetation carbon (Luyssaert et al., 2007), LPJmL5.0 performs well, especially for GPP and NPP, but with a tendency to underestimate vegetation carbon (Fig. 9). The correlation of observed and simulated GPP and NPP values is statistically significant, and the values of vegetation carbon are less correlated owing to the often broad spread of observations within one grid cell (error bars) and general difficulties to exactly reproduce disturbances, mortality, and age class distributions (Schaphoff et al., 2018a). Still, the comparison shows that simulated values are of the right order of magnitude and are also often within the range of observations (error bar crosses 1 : 1 line in Fig. 9). The best correlation with observational data can be found for GPP. We also provide comparisons to eddy flux tower measurements (ORNL DAAC, 2011) in the Supplement. Figures S5–S11 show the modeled versus observed net ecosystem exchange (NEE) rate defined as

$$\begin{array}{}\text{(70)}& \mathrm{NEE}=R_{\text{h}}-\mathrm{NPP},\end{array}$$

where R_h is the heterotrophic respiration. For some sites a time lag (e.g., site Renon–Ritten) between modeled and observed is visible. Because LPJmL5.0 uses the phenology scheme of LPJmL3.5 incorporating the new phenology of LPJmL4 might reduce these deviations (for a comparison, see the Supplement of Schaphoff et al., 2018a). The overall agreement between modeled and observed NEE is satisfying. While LPJmL4 has an averaged Willmott coefficient of agreement (Willmott, 1982) of $\stackrel{\mathrm{‾}}{W}=\mathrm{0.63}$, LPJmL5.0 results in $\stackrel{\mathrm{‾}}{W}=\mathrm{0.58}$. The simulated evapotranspiration fluxes shown in Figs. S12–S20 also agree very well with the observations ($\stackrel{\mathrm{‾}}{W}=\mathrm{0.8}$). The study by Schlesinger and Jasechko (2014), which partitions evapotranspiration into transpiration and evaporation, reports that transpiration accounts for 61 % of global evapotranspiration. In the LPJmL5.0 simulations, transpiration accounts for 59 % on average for the years 2000 to 2009.

4.4.3 Crop yields

The implementation of nitrogen limitation also substantially helps to improve the simulation of global patterns of crop productivity. Regional differences in crop productivity had to be calibrated via a scaling of the maximum leaf area index (LAI_max), the harvest index, and the factor for scaling leaf-level photosynthesis to stand level (α_a), as described in Fader et al. (2010), because the LPJmL3.5 version could only distinguish irrigated from rainfed production but not other inputs such as fertilizers. The difference in crop productivity between low- and high-input regions can now largely be reproduced by the model (compare yellow asterisk for uncalibrated with blue asterisk for calibrated LPJmL3.5 simulations with red circle and cross for LPJmL5 in Figs. S21–S24), which are based on the evaluation procedure as described by Müller et al. (2017). In regions where fertilizer inputs and climate conditions are not the only yield-limiting factors, e.g., in regions with poor pest management, an additional calibration of yield levels could be performed as described by Fader et al. (2010), but is not performed here. The temporal variance of simulated crop yields is often not affected much by accounting for N limitation and sometimes improves or worsens the time series correlation with FAO yield statistics (FAOSTAT, 2018) (Figs. 10 and S25–S27 in the Supplement).

We also use the online tool as supplied by Müller et al. (2017) for comparing the crop yield simulations against the Global Gridded Crop Model Intercomparison (GGCMI) ensemble. Also here, results show that LPJmL5.0 improves with respect to reproducing absolute yield levels across different countries, but there is little effect on the simulated interannual variability of crop yields. As with the calibrated LPJmL3.5 version (Fader et al., 2010), the uncalibrated LPJmL5.0 simulations perform well in comparison to the other GGCMI models. We supply the output of that online model evaluation tool in the Supplement.

5 Discussion and conclusions

The current implementation of nitrogen dynamics into LPJmL3.5, forming LPJmL5.0, introduces a missing feature into a unique modeling framework of the terrestrial biosphere. LPJmL5.0 combines natural vegetation dynamics, the full terrestrial hydrology, and managed grasslands and croplands in one consistent framework with the associated carbon, water, and now also nitrogen pools and fluxes. Owing to parallel model development efforts, LPJmL5.0 does not yet include all model features of the first open source version of LPJmL, LPJmL4 (Schaphoff et al., 2018b, a), such as the updated allocation scheme for managed grasslands (Rolinski et al., 2018) and the updated phenology scheme for natural vegetation (Forkel et al., 2014).

With the implementation of nitrogen dynamics, the model simulations require new inputs, especially on atmospheric deposition, but also on fertilizer applications for which we currently use a static crop- and irrigation-specific data set developed for the harmonization of crop models in the Agricultural Model Intercomparison and Improvement Project (AgMIP) (Elliott et al., 2015). This static fertilizer set also affects the simulation of historic carbon cycle dynamics, as high-input regions such as large parts of Europe and northern America receive current high N inputs in the early 20th century, whereas historic land expansion mostly moves into regions with currently lower input systems. As a consequence, land use change leads to increasing nitrogen limitation and increasing C : N ratios, which may be an artifact from the static fertilizer input data set used.

As historic land use development and fertilizer application are important for simulated current biogeochemical cycles, historic time series of crop-specific fertilizer application would be desirable. Also, global data sets on crop rotations (Kollas et al., 2015), the timing of field operations (Hutchings et al., 2012), and crop residue management, as well as livestock management systems (Rolinski et al., 2018) would be an asset, as the interaction of different cropping systems and natural vegetation is now further increased via the nitrogen cycle.

This first implementation of nitrogen dynamics into LPJmL constitutes an operational modeling framework with many detailed processes resolved explicitly. Specific processes are currently implemented in a simplified manner, even though more detailed approaches are available, such as for biological N fixation (Fisher et al., 2010). As such, some process may have to be revised upon further testing against new reference data. Additional reference data would greatly help to evaluate model performance, which currently is largely constrained to comparisons against other modeling results.

LPJmL5.0 constitutes a unique modeling framework that can now simulate global terrestrial carbon, water, and nitrogen dynamics, consistently accounting for natural vegetation dynamics, agricultural cropland and grassland management, and water management.