A novel representation of biological nitrogen fixation and competitive dynamics between nitrogen-fixing and non-fixing plants in a land model (GFDL LM4.1-BNF)

Representing biological nitrogen fixation (BNF) is an important challenge for coupled carbon (C) and nitrogen (N) land models. Initial representations of BNF in land models applied simplified phenomenological relationships. More recent representations of BNF are mechanistic and include the dynamic response of symbiotic BNF to N limitation of plant growth. However, they generally do not include the competitive dynamics between N-fixing and non-fixing plants, which is a key ecological mechanism that determines ecosystem-scale symbiotic BNF. Furthermore, asymbiotic BNF is generally not included in land models. Here, we present LM4.1-BNF, a novel representation of BNF (asymbiotic and symbiotic) and an updated representation of N cycling in the Geophysical Fluid Dynamics Laboratory Land Model 4.1 (LM4.1). LM4.1-BNF incorporates a mechanistic representation of asymbiotic BNF by soil microbes, a representation of the competitive dynamics between N-fixing and non-fixing plants, and distinct asymbiotic and symbiotic BNF temperature responses derived from corresponding observations. LM4.1-BNF makes reasonable estimations of major carbon (C) and N pools and fluxes and their temporal dynamics, in comparison to the previous version of LM4.1 with N cycling (LM3-SNAP) and to previous representations of BNF in land models generally (phenomenological representations and those without competitive dynamics between N-fixing and non-fixing plants and/or asymbiotic BNF) at a temperate forest site. LM4.1-BNF effectively reproduces asymbiotic BNF rate (13 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in comparison to observations (11 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). LM4.1-BNF effectively reproduces the temporal dynamics of symbiotic BNF rate: LM4.1-BNF simulates a symbiotic BNF pulse in early succession that reaches 73 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 15 years and then declines to ∼0 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 300 years, similarly to observed symbiotic BNF, which reaches 75 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 17 years and then declines to ∼0 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in late successional forests. As such, LM4.1-BNF can be applied to project the dynamic response of vegetation to N limitation of plant growth and the degree to which this will constrain the terrestrial C sink under elevated atmospheric CO₂ concentration and other global change factors.

1 Introduction

The terrestrial carbon (C) sink is controlled by the availability of nitrogen (N) for plant growth (Elser et al., 2007; LeBauer and Treseder, 2008; Wright et al., 2018). Land models are applied to project the terrestrial C sink (Arora et al., 2020) and are progressively incorporating representations of N cycling and N limitation of plant growth (Goll et al., 2017; Lawrence et al., 2019; Medvigy et al., 2009; Smith et al., 2014; Wang et al., 2010; Yang et al., 2009; Zaehle and Friend, 2010). However, the degree to which N limitation of plant growth will constrain the terrestrial C sink under elevated atmospheric CO₂ concentration is unresolved (Terrer et al., 2019), as there is substantial variation between different land models (Wieder et al., 2015b).

The representation of biological N fixation (BNF), the primary natural input of N to terrestrial ecosystems (Fowler et al., 2013; Vitousek et al., 2013), is a key challenge to incorporating N cycling into land models because of its complexity (Davies-Barnard et al., 2020; Meyerholt et al., 2020; Stocker et al., 2016; Thomas et al., 2015; Wieder et al., 2015a). BNF occurs in multiple niches across terrestrial ecosystems: by symbioses between N-fixing bacteria living in root nodules of plants (hereafter, symbiotic BNF) and by a host of other organisms such as soil microbes, bryophytes, and lichens (hereafter, asymbiotic BNF for simplicity although some of these organisms are symbiotic associations; see Reed et al., 2011). Symbiotic BNF and asymbiotic BNF are regulated by a myriad of abiotic and biotic controls, which vary temporally, spatially, and among different niches (Zheng et al., 2019). In particular, symbiotic BNF responds dynamically to N limitation of plant growth: it is up-regulated under N limitation of plant growth and down-regulated under non-N limitation of plant growth (Vitousek et al., 2013). BNF could, as such, be pivotal to overcoming N limitation of plant growth under elevated atmospheric CO₂ concentration (Liang et al., 2016; Terrer et al., 2016, 2018).

Many coupled C–N land models use the empirical relationship of BNF with either net primary production (NPP; Goll et al., 2017) or evapotranspiration (ET; B. Smith et al., 2014; Yang et al., 2009; Zaehle and Friend, 2010) to represent BNF. However, these are simplified phenomenological relationships that are not based on the ecological mechanisms underlying BNF (Cleveland et al., 1999). Furthermore, implementing and comparing a NPP-based and ET-based representation of BNF within a land model (CLM5) resulted in projections of the terrestrial C sink that differed by 50 Pg C in 2100 under the Representative Concentration Pathway 8.5 (RCP8.5; Wieder et al., 2015a). Finally, a recent meta-analysis of BNF found no evidence for the empirical relationship of BNF with either NPP or ET (Davies-Barnard and Friedlingstein, 2020).

Recent coupled C–N land models have simulated symbiotic BNF mechanistically rather than phenomenologically as responding dynamically to N limitation of plant growth. The Geophysical Fluid Dynamics Laboratory (GFDL) Land Model 3 (LM3) can include the Symbiotic Nitrogen Acquisition by Plants (SNAP) model (Sulman et al., 2019), in which plant C allocation to N-fixing bacteria is optimized to maximize plant growth. However, LM3-SNAP and other land models that have implemented a mechanistic representation of symbiotic BNF, such as CLM5 (Lawrence et al., 2019), CABLE (Haverd et al., 2018; Peng et al., 2020; Wang et al., 2010), and E3SM (Zhu et al., 2019), represent a single general plant C pool capable of BNF and cannot represent community dynamics. In observed ecosystems, symbiotic BNF responds dynamically to N limitation of plant growth at both the population scale (via individual-scale regulation of symbiotic BNF rate; Menge et al., 2015) and at the community scale (via competitive dynamics between N-fixing and non-fixing plants; Boring and Swank, 1984; Chapin III et al., 1994; Menge and Hedin, 2009). Under strong N limitation, N-fixing plants up-regulate symbiotic BNF rate and have a competitive advantage over non-fixing plants, but, under weak N limitation, N-fixing plants down-regulate symbiotic BNF rate and are competitively excluded by non-fixing plants because of the high C cost of symbiotic BNF (Gutschick, 1981; Sheffer et al., 2015). As such, the competitive dynamics between N-fixing and non-fixing plants is a key ecological mechanism that could determine ecosystem-scale symbiotic BNF. Finally, the abundance of N-fixing trees is spatially variable (Menge et al., 2019; Staccone et al., 2020), but its representation is not possible in land models that represent a single general plant C pool capable of BNF, although it is necessary to accurately estimate regional symbiotic BNF.

Asymbiotic BNF is generally not included in coupled C–N land models. However, asymbiotic BNF is an important natural input of N to terrestrial ecosystems: in some terrestrial ecosystems, asymbiotic BNF is on par with symbiotic BNF, and asymbiotic BNF has been suggested to account for a substantial proportion of global BNF (Reed et al., 2011). Phenomenological representations of BNF merge asymbiotic and symbiotic BNF, although they are regulated by different controls (Zheng et al., 2019). Mechanistic representations of BNF merge asymbiotic and symbiotic BNF (e.g., LM3-SNAP; Sulman et al., 2019), represent asymbiotic BNF as a constant from averaged observations (e.g., CABLE; Wang and Houlton, 2009), or represent asymbiotic BNF phenomenologically as a function of ET (e.g., CLM5; Lawrence et al., 2019). Importantly, although asymbiotic and symbiotic BNF exhibit different temperature responses (Bytnerowicz et al., 2021), the symbiotic BNF temperature response is, when included, derived primarily from asymbiotic BNF observations (Houlton et al., 2008), and the asymbiotic BNF temperature response is omitted.

Here, we present LM4.1-BNF, a novel representation of BNF and an updated representation of N cycling in the GFDL land model 4.1 (LM4.1; Shevliakova et al., 2021). LM4.1 includes height-structured competition for light and water between plant cohorts using the perfect plasticity approximation (Martinez Cano et al., 2020; Purves et al., 2008; Strigul et al., 2008; Weng et al., 2015). LM4.1-BNF builds on the framework of LM4.1, including competition for light, water, and N between plant cohorts that associate with N-fixing bacteria and non-fixer plant cohorts. LM4.1-BNF introduces several improvements to the representation of N cycling in LM3-SNAP by incorporating novel representations of the following ecological mechanisms.

1. Symbiotic BNF and competitive dynamics between N-fixing and non-fixing plants. Plant cohorts with a N-fixer vegetation type conduct symbiotic BNF and compete with plant cohorts with a non-fixer vegetation type.

2. Asymbiotic BNF. Soil microbes conduct asymbiotic BNF, as well as decomposition, nitrification, and denitrification.

3. BNF temperature response. Asymbiotic and symbiotic BNF have different temperature responses derived from asymbiotic BNF observations (Houlton et al., 2008) and symbiotic BNF observations (Bytnerowicz et al., 2021) respectively.

4. N limitation. N limitation is determined by current stored non-structural N relative to the demand for non-structural N. N limitation increases active root uptake of inorganic N and decreases root N exudation following observations (Canarini et al., 2019; Nacry et al., 2013).

5. Dynamic plant C allocation to growth and N uptake. N limitation decreases the growth of leaves, sapwood, and seeds, proportionally increasing the growth of fine roots following observations (Poorter et al., 2012). N limitation stimulates C allocation to N uptake (including symbiotic BNF) relative to growth. C limitation, which is determined by current stored non-structural C relative to the demand for non-structural C, stimulates C allocation to growth relative to N uptake. Thereby, plant C allocation is optimized to maximize growth following observations (Rastetter and Shaver, 1992).

We focus our analysis on temperate forests which are generally N-limited (Elser et al., 2007; LeBauer and Treseder, 2008). We parameterize a N-fixer vegetation type based on Robinia pseudoacacia (black locust), which is the most abundant N-fixing tree species in the coterminous United States, accounting for 64 % of tree-associated BNF in the coterminous United States (Staccone et al., 2020). We compare Robinia to a non-fixer vegetation type based on Acer rubrum (red maple), which is the most abundant non-fixing tree species in the north region of the coterminous United States (Oswalt et al., 2019). We evaluate LM4.1-BNF at Coweeta Hydrologic Laboratory in North Carolina, United States, which has observations on symbiotic BNF by Robinia (Boring and Swank, 1984).

We conduct three analyses to assess the performance of LM4.1-BNF in estimating major C and N pools and fluxes in comparison to previous representations of BNF in land models generally. In the first analysis, we compare mechanistic and phenomenological representations of BNF. We compare LM4.1-BNF (with BNF represented mechanistically as described above) to LM4.1-BNF with BNF represented as a function of NPP and to LM4.1-BNF with BNF represented as a function of ET. In the second analysis, we examine the role of competitive dynamics between N-fixing and non-fixing plants. We compare LM4.1-BNF simulations with both Robinia and Acer to LM4.1-BNF simulations with only Acer and LM4.1-BNF simulations with only Acer that can associate with N-fixing bacteria, which are representative of land models that represent a single general plant C pool capable of BNF and cannot represent community dynamics. In the third analysis, we examine the role of asymbiotic BNF. We compare LM4.1-BNF simulations with asymbiotic BNF to LM4.1-BNF simulations without asymbiotic BNF, which is representative of land models that do not include asymbiotic BNF.

2 Model description

2.1 Overview of a land tile and vegetation types

We provide an overview of LM4.1-BNF with a focus on the novel elements relative to LM4.1 (Shevliakova et al., 2021) and LM3-SNAP (Sulman et al., 2019). A complete description of LM4.1-BNF is available in Appendix A. Note that LM4.1 can be coupled with the GFDL atmosphere model to serve as a base for the GFDL climate and Earth system models (Zhao et al., 2018a, b).

LM4.1-BNF consists of a grid, in which grid cells are approximately 100 km by 100 km. LM4.1-BNF represents the heterogeneity of the land surface as a mosaic of land tiles within a grid cell. Each land tile represents a fraction of the grid cell area and does not have an associated location within the grid cell. A land tile may represent natural vegetation at a given stage of recovery post-disturbance, urban area, pastureland, rangeland, or cropland. Land tiles are created dynamically due to a disturbance, such as human land use, fire, or natural mortality of vegetation.

A land tile contains multiple plant cohorts that compete for light and water following the perfect plasticity approximation (Martinez Cano et al., 2020; Purves et al., 2008; Strigul et al., 2008; Weng et al., 2015) and compete for N (presented below). Plant cohorts consist of identical individual trees belonging to a vegetation type that occupy a given canopy layer and that have a spatial density (determined by recruitment and mortality). A vegetation type can be associated with exclusively arbuscular mycorrhizae (AM), exclusively ectomycorrhizae (EM), both AM and N-fixing bacteria, or both EM and N-fixing bacteria. A land tile can contain multiple plant cohorts of the same or of different vegetation types. As such, there is intraspecific competition (among plant cohorts of the same vegetation type within a tile) and interspecific competition (among plant cohorts of different vegetation types within a tile). Growth is based on allometric equations (Eqs. A40–A42) and is modulated by N availability. Recruitment and mortality follow Weng et al. (2015) and Martinez Cano et al. (2020) and are not directly influenced by N availability but are indirectly influenced by N availability via its effect on growth.

Table 1 Key parameter differences between vegetation types. See Table D1 for remaining vegetation-type-specific parameters.

There are six plant tissue C and N pools: leaf, fine root, sapwood, heartwood, seed, and non-structural C or N. The C:N ratios of the leaf, fine root, sapwood, heartwood, and seed tissue pools are fixed (for a given vegetation type) (Tables 1, D1, and D2). There are three soil organic C and N pools (labile plant-derived, labile microbe-derived, and recalcitrant) and two soil inorganic N pools (ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$)) in each soil layer. There are 20 soil layers of varying thickness to a total depth of 10 m. Soil C and N are transferred between soil layers via leaching. The soil C:N ratio is not fixed. Figure 1 displays a diagram of key C and N pools and fluxes.

Figure 1 Diagram of key C and N pools (boxes) and fluxes (arrows). C pools and fluxes are indicated in blue. N pools and fluxes are indicated in orange. NSC represents non-structural C and NSN represents non-structural N. Orange and blue boxes have a fixed C:N ratio. Plant turnover, symbiont turnover, and soil C and N pools and fluxes are not displayed for visual clarity. The black dashed box represents the rhizosphere. (a) Vegetation type that associates with arbuscular mycorrhizae (AM). (b) Vegetation type that associates with AM and N-fixing bacteria. (c) Vegetation type that associates with ectomycorrhizae (EM). (d) Vegetation type that associates with EM and N-fixing bacteria.

We define a N-fixer vegetation type with a parameterization based on Robinia pseudoacacia and a non-fixer vegetation type with a parameterization based on Acer rubrum. Both Acer and Robinia associate with AM. We used the US Forest Inventory and Analysis (FIA) database (US Forest Service, 2020a), the US FIA Forest Health Monitoring database (US Forest Service, 2020b), and the Biomass and Allometry Database (BAAD; Falster et al., 2015) to parameterize the allometries of these vegetation types (Eqs. A40–A42; Appendices B and C). The vegetation types also differed in other key traits (Table 1). In particular, the C:N ratio of leaves differed between vegetation types that associated with AM, EM, and N-fixing bacteria (Adams et al., 2016; Averill et al., 2019). See Table D1 for all vegetation-type-specific parameters. Other model parameters are from Weng et al. (2015) or Sulman et al. (2019) or are derived from published observations (Appendix C). Some parameters were not well constrained by available observations and were tuned to fit to observed patterns of C and N cycling in temperate forests (Appendix C). See Table D2 for general parameters.

2.2 Symbiotic BNF and N uptake by roots, AM, and EM

All vegetation types take up inorganic N via passive and active root uptake. Passive root uptake of inorganic N follows LM3-SNAP (Eq. A1). Active root uptake of inorganic N follows LM3-SNAP but is modified to increase with N stress following observations (Nacry et al., 2013) (described in further detail below; Eq. 16). AM takes up inorganic N following LM3-SNAP (Eq. A4). EM decomposes and takes up organic C and N following LM3-SNAP but is modified to additionally take up inorganic N following observations (Phillips et al., 2013) (Eq. A8).

The symbiotic BNF rate by N-fixing bacteria (N_Nfix; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(1)}& {\mathrm{N}}_{\text{Nfix}}=r_{\text{Nfix}}{B}_{\text{Nfix}}{f}_{\mathrm{s}}\left(T\right),\end{array}$$

where r_Nfix is a rate constant, B_Nfix is the biomass C of the nodule (includes both plant and N-fixing bacteria tissue) [kg N indiv⁻¹], and f_s(T) is the soil temperature dependence function. For Robinia,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{f}_{\mathrm{s}}\left(T\right)=& max[\mathrm{0.0},\left({\displaystyle \frac{\mathrm{45.67}-\left(T-\mathrm{273.15}\right)}{\mathrm{45.67}-\mathrm{31.89}}}\right)\\ & {\left({\displaystyle \frac{\left(T-\mathrm{273.15}\right)-\mathrm{1.43}}{\mathrm{31.89}-\mathrm{1.43}}}\right)}^{\frac{\mathrm{31.89}-\mathrm{1.43}}{\mathrm{45.67}-\mathrm{31.89}}}],\end{array}\end{array}$$

where T is the average soil temperature across soil layers [K]. This function reaches its maximum at 31.9 ^∘C (Fig. 2). It is a modified beta distribution function (Yan and Hunt, 1999) and is derived from Bytnerowicz et al. (2021).

Figure 2 Soil temperature dependence function of asymbiotic BNF and symbiotic BNF by Robinia. Asymbiotic BNF reaches its maximum at 24.4 ^∘C and is derived from Houlton et al. (2008). Symbiotic BNF by Robinia reaches its maximum at 31.9 ^∘C and is derived from Bytnerowicz et al. (2021).

Respiration associated with symbiotic BNF is cost_NfixN_Nfix, where cost_Nfix is the C cost of symbiotic BNF per unit N (Eq. A52).

Note that the description of EM is included although both Acer and Robinia associate with AM.

2.3 Asymbiotic BNF

Soil microbes are represented as a single C pool that conducts decomposition, nitrification, denitrification, and asymbiotic BNF. The rates of C and N decomposition, rates of C and N decomposition during denitrification, rates of change of biomass C and N, and maintenance respiration rate of soil microbes follow LM3-SNAP (Eqs. A11–A19). The N surplus or deficit of soil microbes and C and N growth rates of soil microbes are modified to include asymbiotic BNF (Eqs. A22–A24).

The asymbiotic BNF rate of soil microbes in soil layer k (N_{Nfix asymb}(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(3)}& {\mathrm{N}}_{\text{Nfix asymb}}\left(k\right)=r_{\text{Nfix asymb}}{\mathrm{C}}_{\mathrm{M}}\left(k\right)f_a\left(T\left(k\right)\right),\end{array}$$

where r_{Nfix asymb} is a rate constant, C_M(k) is the biomass C of soil microbes in soil layer k [kg C m⁻²], and f_a(T(k)) is the soil temperature dependence function.

$$\begin{array}{}\text{(4)}& f_a\left(T\left(k\right)\right)=e^{-\mathrm{2.6}+\mathrm{0.21}\left(T\left(k\right)-\mathrm{273.15}\right)\left(\mathrm{1}-\frac{\mathrm{0.5}\left(T\left(k\right)-\mathrm{273.15}\right)}{\mathrm{24.4}}\right)}\end{array}$$

This function reaches its maximum at 24.4 ^∘C (Fig. 2). It is a modified normal distribution function and is derived from the observations compiled by Houlton et al. (2008) with the study of symbiotic BNF removed (Schomberg and Weaver, 1992) and is normalized to a maximum of 1.

2.4 N limitation, plant C allocation to growth, and plant C allocation to rhizosphere priming

The non-structural C pool (NSC; [kg C indiv⁻¹]) gains C from photosynthesis. NSC loses C to respiration and C allocation to growth, symbionts, and root C exudation. The rate of change of NSC ($\frac{\text{dNSC}}{\mathrm{d}t}$; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle \frac{\text{dNSC}}{\mathrm{d}t}}=& {\displaystyle }P-R-\left(G_{\mathrm{C},\mathrm{l}}+G_{\mathrm{C},\mathrm{r}}+G_{\mathrm{C},\text{sw}}+G_{\mathrm{C},\text{seed}}\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }-{\mathrm{C}}_{\text{alloc}}-L_{\mathrm{C},\text{exudate}},\end{array}$$

where P is the photosynthesis rate [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], R is the respiration rate (maintenance and growth) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,l is the growth rate of the leaf C pool (C_l; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,r is the growth rate of the fine-root C pool (C_r; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], C_sw is the growth rate of the sapwood C pool (C_sw; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,seed is growth rate of the seed C pool (C_seed; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], C_alloc is the rate of C allocation to symbionts (Eqs. A43–A46), and L_C,exudate is the rate of root C exudation (Eq. A38). Note that sapwood is converted to heartwood following Martinez Cano et al. (2020).

The non-structural N pool (NSN; [kg N indiv⁻¹]) gains N from N uptake via roots and symbionts. NSN loses N to N allocation to growth, symbionts, and root N exudation. The rate of change of NSN ($\frac{\text{dNSN}}{\mathrm{d}t}$; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle \frac{\text{dNSN}}{\mathrm{d}t}}=& {\displaystyle }U-\left({\displaystyle \frac{G_{\mathrm{C},\mathrm{l}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{l}}}}+{\displaystyle \frac{G_{\mathrm{C},\mathrm{r}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{r}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{sw}}}{{\mathrm{C}:\mathrm{N}}_{\text{sw}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{seed}}}{{\mathrm{C}:\mathrm{N}}_{\text{seed}}}}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }-{\mathrm{N}}_{\text{alloc}}-L_{\mathrm{N},\text{exudate}},\end{array}$$

where U is the N uptake rate via roots and symbionts (Eq. A55) [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], C:N_l is the fixed C:N ratio of leaves, C:N_r is the fixed C:N ratio of fine roots, C:N_sw is the fixed C:N ratio of sapwood, C:N_seed is the fixed C:N ratio of seeds, N_alloc is the rate of N allocation to symbionts (Eq. A47), and L_N,exudate is the rate of root N exudation (Eq. A39).

Non-N-limited growth is calculated according to Weng et al. (2015) and Martinez Cano et al. (2020). The total allocation of NSC to growth is determined by the target NSC (NSC_target; [kg C indiv⁻¹]) and minimizes the deviation between NSC and NSC_target. NSC_target is a multiple of the target C_l (C_l,target; [kg C indiv⁻¹]), which reflects the ability of a plant to refoliate after defoliation (Hoch et al., 2003; Richardson et al., 2013), and is calculated as

$$\begin{array}{}\text{(7)}& {\text{NSC}}_{\text{target}}=q{\mathrm{C}}_{\mathrm{l},\text{target}},\end{array}$$

where q is a proportionality constant. The allocation of NSC to the growth of each tissue depends on the total allocation of NSC to growth and the target C pool of each tissue and minimizes the deviation between the C pool of each tissue and the target C pool of each tissue. The target C pool of each tissue is dynamic and is determined by allometry (Eqs. A40–A42), canopy position, and phenology.

In LM4.1-BNF, G_C,l, G_C,r, G_C,sw, and G_C,seed are adjusted to include N limitation and are calculated as

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\mathrm{l}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\mathrm{l}},\text{(9)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\mathrm{r}}={\mathrm{\Delta }}_{\mathrm{r}},\text{(10)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\text{sw}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\text{sw}},\text{(11)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\text{seed}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\text{seed}},\end{array}$$

where N_stress is N stress [unitless] and Δ_l, Δ_r, Δ_sw, and Δ_seed are the non-N-limited growth rates of C_l, C_r, C_sw, and C_seed respectively [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$] following Weng et al. (2015) and Martinez Cano et al. (2020). Because plants increase C allocation to fine roots relative to other tissues when N-limited (Poorter et al., 2012), G_C,r is not adjusted to include N limitation.

In LM4.1-BNF, N_stress is the relative difference between NSN and NSN_target and is calculated as

$$\begin{array}{}\text{(12)}& {\mathrm{N}}_{\text{stress}}=max\left[\mathrm{0},{\displaystyle \frac{{\text{NSN}}_{\text{target}}-\text{NSN}}{{\text{NSN}}_{\text{target}}}}\right],\end{array}$$

where NSN_target is the target NSN [kg N indiv⁻¹]. N_stress is smoothed with a low-pass filter over 30 d to reflect the persisting influence of N stress (Mooney et al., 1991). NSN_target is calculated as

$$\begin{array}{}\text{(13)}& {\text{NSN}}_{\text{target}}={\displaystyle \frac{{\text{NSC}}_{\text{target}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{l}}}}.\end{array}$$

This is similar to LM3-SNAP, which compared the target leaf and root N pools to NSN, but is modified to reflect the treatment of NSC_target in LM4.1 by including the target sapwood and seed N pools.

Plant turnover decreases C_l, C_r, and C_sw and from N_l, N_r, and N_sw at a constant tissue-specific rate. A fraction of the turnover of C_l and N_l is retranslocated into NSC and NSN respectively.

Under N limitation, plants increase root C exudation to stimulate N mineralization in the rhizosphere (rhizosphere priming; Cheng et al., 2014; Finzi et al., 2015). L_C,exudate increases with N_stress and is calculated as

$$\begin{array}{}\text{(14)}& L_{\mathrm{C},\text{exudate}}=r_{\text{leakage},\mathrm{C}}\text{NSC}{\mathrm{N}}_{\mathrm{stress}},\end{array}$$

where r_leakage,C is a rate constant.

Under N limitation, plants decrease root N exudation (Canarini et al., 2019). L_N,exudate decreases with N_stress and is calculated as

$$\begin{array}{}\text{(15)}& L_{\mathrm{N},\text{exudate}}=r_{\text{leakage},\mathrm{N}}\text{NSN}\left(\mathrm{1}-{\mathrm{N}}_{\text{stress}}\right),\end{array}$$

where r_leakage,N is a rate constant.

Under N limitation, plants increase active root uptake of inorganic N (Nacry et al., 2013). The rate of active root uptake of inorganic N in soil layer k (N_active(k); [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) increases with N_stress and is calculated as

$$\begin{array}{}\text{(16)}& {\mathrm{N}}_{\text{active}}\left(k\right)=f\left({\mathrm{NO}}_{\mathrm{3}}\left(k\right),{\mathrm{NH}}_{\mathrm{4}}\left(k\right),f_{\text{rhiz}}\left(k\right)\right){\mathrm{N}}_{\text{stress}},\end{array}$$

where $f\left({\mathrm{NO}}_{\mathrm{3}}\right(k),{\mathrm{NH}}_{\mathrm{4}}(k),f_{\text{rhiz}}(k\left)\right)$ is a function of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in soil layer k (NH₄(k); [kg N m⁻²]), the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool in soil layer k (NO₃(k); [kg N m⁻²]), and the rhizosphere volume fraction of soil layer k (f_rhiz(k); [m³ m⁻³]) and is given in Eq. (A3).

2.5 Plant C allocation to symbionts (AM, EM, and N-fixing bacteria)

The rate of C allocation to AM (C_alloc,AM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(17)}& {\mathrm{C}}_{\text{alloc},\text{AM}}=f_{\text{alloc},\text{AM}}\text{NSC},\end{array}$$

where f_alloc,AM is the fraction of NSC allocated to AM per unit time. C_alloc,AM is not related to N_stress because, although AM increases N uptake, AM is maintained by the plant primarily for phosphorus uptake (Smith and Smith, 2011).

The rate of C allocation to EM (C_alloc,EM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(18)}& {\mathrm{C}}_{\text{alloc},\text{EM}}=f_{\text{alloc},\text{EM}}\text{NSC}{\mathrm{N}}_{\text{stress}},\end{array}$$

where f_alloc,EM is the maximum fraction of NSC allocated to EM per unit time. C_alloc,EM is a function of N_stress because biomass C of EM increases with N limitation (Phillips et al., 2013).

Plants that associate with N-fixing bacteria can regulate symbiotic BNF to different extents, termed their BNF strategy (Menge et al., 2015). For plants with a perfectly facultative BNF strategy, symbiotic BNF increases with N limitation. For plants with an incomplete BNF strategy, symbiotic BNF increases with N limitation but is maintained at a minimum. For plants with an obligate BNF strategy, symbiotic BNF is constant. For plants with either a facultative or an incomplete BNF strategy, the rate of C allocation by the plant to N-fixing bacteria (C_alloc,Nfix; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(19)}& {\mathrm{C}}_{\text{alloc},\text{Nfix}}=max\left[f_{\text{alloc},\text{Nfix}}\text{NSC}{\mathrm{N}}_{\text{stress}},f_{\text{alloc},\text{Nfix},\text{min}}\text{NSC}\right],\end{array}$$

where f_alloc,Nfix is the fraction of NSC allocated to N-fixing bacteria per unit time, and $f_{\text{alloc},\text{Nfix},\text{min}}$ is the minimum fraction of NSC allocated to N-fixing bacteria per unit time. For plants with a perfectly facultative BNF strategy $f_{\text{alloc},\text{Nfix},\text{min}}=\mathrm{0}$, and for plants with an incomplete down-regulator BNF strategy $f_{\text{alloc},\text{Nfix},\text{min}}>\mathrm{0}$. Robinia has an incomplete down-regulator BNF strategy.

For plants with an obligate BNF strategy, C_alloc,Nfix is

$$\begin{array}{}\text{(20)}& {\mathrm{C}}_{\text{alloc},\text{Nfix}}=f_{\text{alloc},\text{Nfix}}\text{NSC}.\end{array}$$

Additionally, plants allocate a small quantity of N to symbionts such that symbiont growth can be initiated. The rate of N allocation by the plant to symbionts (N_alloc,j; j=AM, EM, N-fixing bacteria; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(21)}& {\mathrm{N}}_{\text{alloc},j}={\displaystyle \frac{{\mathrm{C}}_{\text{alloc},j}}{{\mathrm{C}:\mathrm{N}}_{\text{alloc}}}},\end{array}$$

where C:N_alloc is the C:N ratio of C and N allocated to symbionts by the plant.

Plant C allocation to symbionts increases biomass C of symbionts, which increases plant N uptake via symbionts (Appendices A1 and A6).

2.6 Dynamic plant C allocation to growth and N uptake

The order of plant C allocation to growth, symbionts, and rhizosphere priming is determined by C limitation relative to N limitation (Cheng et al., 2014; Finzi et al., 2015; Poorter et al., 2012; Treseder, 2004; Zheng et al., 2019). If a plant is more C-limited than N-limited, $\text{NSC}<\text{NSN}\cdot {\mathrm{C}:\mathrm{N}}_{\text{leaf}}$, the plant allocates C to growth, then to N-fixing bacteria (if associated) and EM (if associated), and then to rhizosphere priming, and finally to AM. If a plant is more N-limited than C-limited, $\text{NSC}>\text{NSN}\cdot {\mathrm{C}:\mathrm{N}}_{\text{leaf}}$, the plant allocates C to N-fixing bacteria (if associated) and EM (if associated), then to rhizosphere priming, then to growth, and finally to AM.

2.7 Soil N₂O and NO emissions

Soil N₂O and NO emissions occur during nitrification (aerobic oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ with oxygen as an electron acceptor, which produces N₂O and NO as by-products) and denitrification (anaerobic oxidation of organic C with ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as an electron acceptor, which produces N₂O as a by-product).

Following LM3V-N (Huang and Gerber, 2015), soil N₂O emission rate in soil layer k (N₂O(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\left(k\right)=& {\displaystyle }\mathrm{0.004}\text{nit}\left(k\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{1}/[\mathrm{1}+max\left(\mathrm{3.52},\mathrm{22}\exp \left(-\mathrm{0.8}{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\text{HR}\left(k\right)}}\right)\right)\\ \text{(22)}& {\displaystyle }& {\displaystyle }\cdot max\left(\mathrm{0.1},\mathrm{0.015}{\displaystyle \frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}}-\mathrm{0.32}\right)\text{denit}\left(k\right)],\end{array}$$

where nit(k) is nitrification rate in soil layer k [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] (Eq. A56), HR(k) is heterotrophic respiration in soil layer k [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$], θ(k) is volumetric soil water content of soil layer k [m³ m⁻³], θ_sat is saturation volumetric soil water content, and denit(k) is denitrification rate in soil layer k [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] (Eq. A57).

Following LM3V-N (Huang and Gerber, 2015), soil NO emission rate in soil layer k (NO(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{NO}\left(k\right)=\mathrm{0.004}\left(\mathrm{15.2}+{\displaystyle \frac{\mathrm{35.5}{\tan }^{-\mathrm{1}}\left(\mathrm{0.68}\mathit{\pi }\left(\mathrm{2.09}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{4}/\mathrm{3}}-\mathrm{1.68}\right)\right)}{\mathit{\pi }}}\right)\\ \text{(23)}& {\displaystyle }& {\displaystyle }\cdot \text{nit}\left(k\right).\end{array}$$

2.8 Phenomenological representations of BNF (LM4.1-BNF_NPP and LM4.1-BNF_ET)

To compare mechanistic and phenomenological representations of BNF in our first analysis, we developed LM4.1-BNF_NPP and LM4.1-BNF_ET which have phenomenological representations of BNF. In LM4.1-BNF_NPP, BNF is represented as a function of NPP in LM4.1-BNF. BNF rate (BNF_NPP; [kg N m⁻²]) is

$$\begin{array}{}\text{(24)}& {\text{BNF}}_{\text{NPP}}=a_{\text{NPP}}\left(\mathrm{1}-e^{b_{\text{NPP}}\text{NPP}}\right),\end{array}$$

where NPP is net primary production [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$], and a_NPP and b_NPP are constants from Meyerholt et al. (2016).

In LM4.1-BNF_ET, BNF is represented as a function of ET in LM4.1-BNF. BNF rate (BNF_ET; [kg N m⁻²]) is

$$\begin{array}{}\text{(25)}& {\text{BNF}}_{\text{ET}}=a_{\text{ET}}\text{ET}+b_{\text{ET}},\end{array}$$

where ET is evapotranspiration [mm yr⁻¹], and a_ET and b_ET are constants from Meyerholt et al. (2016).

BNF_NPP and BNF_ET enter NH₄(k) (distributed across all soil layers proportional to thickness). In LM4.1-BNF_NPP and LM4.1-BNF_ET, growth and turnover of symbionts and plant C allocation to symbionts do not occur, and asymbiotic BNF does not occur. All other components of C and N cycling in LM4.1-BNF_NPP and LM4.1-BNF_ET are the same as LM4.1-BNF.

3 Numerical experiments and evaluation description

3.1 Numerical experiments description

We ran numerical experiments for the grid cell containing Coweeta Hydrological Laboratory (CHL) in North Carolina, United States (35.05^∘ N, 83.45^∘ W), which is part of the Long-Term Ecological Research Network and has observations on symbiotic BNF by Robinia (Boring and Swank, 1984).

We ran the LM4.1-BNF spin-up for 1000 years at pre-industrial atmospheric CO₂ concentration (284.26 ppm) to allow the soil C and N pools to reach an approximate steady state. Then, we initialized LM4.1-BNF, LM4.1-BNF_NPP, and LM4.1-BNF_ET numerical experiments with seedlings (removed vegetation C and N pools from the LM4.1-BNF spin-up) and the LM4.1-BNF spin-up soil C and N pools. We ran numerical experiments for another 300 years at current atmospheric CO₂ concentration (324.53 ppm). See Table D3 for a summary of atmospheric CO₂ concentration (Dlugokencky and Tans, 2020; Meinshausen et al., 2017), meteorological forcings (Sheffield et al., 2006), and N deposition rates (Dentener, 2006) used in the spin-up and numerical experiments. We initialized the LM4.1-BNF spin-up with Acer seedlings, and we initialized LM4.1-BNF numerical experiments with both Acer and Robinia seedlings.

Table 2 Description of numerical experiments. In the first analysis, we compare mechanistic and phenomenological representations of BNF. In the second analysis, we examine the role of competitive dynamics between N-fixing and non-fixing plants. In the third analysis, we examine the role of asymbiotic BNF. Note that the same six numerical experiments were examined in both the second and third analyses.

To compare mechanistic and phenomenological representations of BNF in our first analysis (Table 2), we initialized LM4.1-BNF_NPP and LM4.1-BNF_ET numerical experiments with only Acer seedlings. To examine the role of competitive dynamics between N-fixing and non-fixing plants in our second analysis (Table 2), we initialized LM4.1-BNF numerical experiments with both Acer and Robinia seedlings, only Acer seedlings, or only Acer seedlings that can associate with N-fixing bacteria (N-fixer Acer). To examine the role of asymbiotic BNF in our third analysis (Table 2), we initialized LM4.1-BNF numerical experiments with both Acer and Robinia seedlings, only Acer seedlings, or only N-fixer Acer seedlings. LM4.1-BNF, LM4.1-BNF_NPP, and LM4.1-BNF_ET simulations are initialized such that all plant cohorts have the same height (0.5 m; Table D4), and dbh is determined from height by allometry (Eq. A41; Table D4). See Table D4 for a summary of initial density, height, and dbh of seedlings in numerical experiments.

We ran the LM3-SNAP spin-up for 1000 years at pre-industrial atmospheric CO₂ concentration to allow the soil C and N pools to reach an approximate steady state. Then, we initialized LM3-SNAP numerical experiments with seedlings and the LM3-SNAP spin-up soil C and N pools. We ran numerical experiments for another 300 years at current atmospheric CO₂ concentration. See Table D3 for a summary of atmospheric CO₂ concentration, meteorological forcings, and N deposition rates used in the spin-up and numerical experiments. We initialized the LM3-SNAP spin-up and the LM3-SNAP numerical experiments with a temperate deciduous vegetation type.

3.2 Evaluation description

To evaluate the performance of LM4.1-BNF in representing symbiotic BNF, we compared symbiotic BNF rate from the numerical experiments of LM4.1-BNF and LM3-SNAP to site observations from CHL (Fig. 3a; Boring and Swank, 1984). To evaluate the performance of LM4.1-BNF in representing asymbiotic BNF, we compared asymbiotic BNF rate from the numerical experiments of LM4.1-BNF to site observations from CHL (Fig. 3b; Todd et al., 1978). To evaluate the performance of LM4.1-BNF in representing the competitive dynamics between N-fixing and non-fixing plants, we compared the basal area fraction of each vegetation type (Acer and Robinia) over time from the numerical experiments of LM4.1-BNF to US FIA database tree data (Fig. 4; US Forest Service, 2020a).

Figure 3 Simulated BNF rate by LM4.1-BNF and LM3-SNAP. (a) Simulated asymbiotic BNF rate compared to CHL site data. Simulated data are averaged over the last 100 years of the 300 years of simulation to reflect the site data which are from mature forests. Error bars indicate 2 standard deviations. (b) Simulated symbiotic BNF rate over time compared to CHL site data for a 4, 17, and 38 years and mature forest (plotted at 300 years).

Figure 4 Simulated relative basal area of Acer and Robinia over time compared to FIA data (in North Carolina). Simulated data are trees with dbh>12.7 cm to reflect the dbh range of FIA data. FIA data of all non-fixing trees are aggregated to represent Acer. Each point represents an FIA plot. See Fig. D9 for absolute basal area.

Figure 5 Simulated total plant biomass C over time by LM4.1-BNF and LM3-SNAP compared to FIA data (in North Carolina). Each point represents an FIA plot.

Figure 6 Simulated soil C and N pools, soil N fluxes, and soil N loss rates by LM4.1-BNF and LM3-SNAP. (a) Simulated total soil C (depth 0–10 cm) compared to FIA data (in North Carolina) and CHL site data. (b) Simulated total soil N (depth 0–10 cm) compared to FIA data (in North Carolina) and CHL site data. (c) Simulated soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (depth 0–10 cm) compared to CHL site data. (d) Simulated N mineralization rate and net nitrification rate (depth 0–10 cm) compared to CHL site data. (e) Simulated N₂O and NO emission rates compared to a meta-analysis estimate for temperate forests and simulated dissolved organic N (DON), ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate compared to CHL site data. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the site data which are from mature forests. Error bars indicate 2 standard deviations. NA indicates that LM3-SNAP cannot estimate N₂O or NO emissions.

To evaluate the performance of LM4.1-BNF in representing ecosystem C and N pools and fluxes, we (1) compared total plant biomass C from the numerical experiments of LM4.1-BNF and LM3-SNAP to US FIA database tree data (Fig. 5; US Forest Service, 2020a); (2) compared soil C and N pools and fluxes from the numerical experiments of LM4.1-BNF and LM3-SNAP to US FIA database soil data (Fig. 6; US Forest Service, 2020a), site observations from CHL (Fig. 6; Knoepp, 2009a, b, 2018; Swank and Waide, 1988), and a meta-analysis of temperate forests (Fig. 6; Stehfest and Bouwman, 2006); (3) compared ecosystem C fluxes from the numerical experiments of LM4.1-BNF and LM3-SNAP to eddy covariance observations from CHL at the hourly timescale (Fig. 7; Oishi, 2020); and (4) compared ecosystem C fluxes from the numerical experiments of LM4.1-BNF and LM3-SNAP to a meta-analysis of temperate forests (Fig. 8; Anderson-Teixeira et al., 2018). We also compared dbh growth rates and the dbh distribution from the numerical experiments of LM4.1-BNF to US FIA database tree data (Figs. D1 and D2; US Forest Service, 2020a). See Appendix B for data availability and processing. See Table D5 for a summary of the validated variables and data sources.

Figure 7 Simulated net ecosystem production (NEP) by LM4.1-BNF and LM3-SNAP compared to net ecosystem exchange (NEE) CHL site data at the hourly timescale. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the data which are from mature forests.

4 Evaluation: comparison of LM4.1-BNF to observations and to LM3-SNAP

4.1 Evaluation results

Here we describe the evaluation of LM4.1-BNF, in which we compare LM4.1-BNF simulations to observations and to LM3-SNAP simulations (which represents a single general plant C pool capable of BNF and cannot represent community dynamics (i.e., competitive dynamics between N-fixing and non-fixing plants) and does not represent asymbiotic BNF). LM4.1-BNF captures observed symbiotic BNF, asymbiotic BNF, successional dynamics, and the major pools and fluxes of C and N.

LM4.1-BNF effectively reproduces asymbiotic BNF rate (mean 13 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the final 100 years), which is not represented in LM3-SNAP (Fig. 3a), in comparison to observations from CHL (11 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). LM4.1-BNF effectively reproduces the temporal dynamics of symbiotic BNF rate: LM4.1-BNF simulates a symbiotic BNF rate pulse in early succession that reaches 73 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 15 years and then declines to ∼0 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 300 years (Fig. 3b). LM3-SNAP simulates high symbiotic BNF rate in late succession (mean 8 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the final 100 years; Fig. 3b). Observations from CHL suggest that symbiotic BNF which reaches 75 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ at 17 years then declines to ∼0 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in late successional forests. Note that asymbiotic BNF is directly controlled by soil microbe biomass C (Eq. 3) and peripherally controlled by soil microbe biomass N, total soil C, and total soil N (Fig. D3), and symbiotic BNF is directly controlled by nodule biomass C (Eq. 1) and peripherally controlled by non-structural C and N stress (Fig. D4).

LM4.1-BNF effectively reproduces the successional dynamics of Robinia and Acer (Fig. 4): Robinia is competitively excluded by Acer at approximately the same timescale as observations (∼5 % basal area fraction at 150 years).

At the ecosystem scale, LM4.1-BNF effectively reproduces the temporal dynamics of total plant biomass C in comparison to observations (mean 173 vs. 160 Mg C ha⁻¹ over the final 100 years; Fig. 5), whereas LM3-SNAP substantially underestimates total plant biomass C in comparison to observations (mean 44 vs. 160 Mg C ha⁻¹ over the final 100 years; Fig. 5). LM4.1-BNF makes reasonable estimates for soil C and N pools and fluxes, which are comparable to those of LM3-SNAP. LM4.1-BNF underestimates total soil C and N (mean 13.4 vs. 33.0 Mg C ha⁻¹ and 0.2 vs. 1.8 Mg N ha⁻¹ respectively; Fig. 6a and b) and overestimates soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (mean 4.7 vs. 1.0 and 1.6 vs. 0.1 kg N ha⁻¹ respectively; Fig. 6c) in comparison to observations. LM4.1-BNF underestimates N mineralization rate and net nitrification rate in comparison to observations (mean 14.9 vs. 39.5 and 7.2 vs. 12.8 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ respectively; Fig. 6d). LM4.1-BNF overestimates N₂O emission rate (mean 5.1 vs. 0.9 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and underestimates NO emission rate (mean 0.2 vs. 0.8 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 6e) in comparison to observations. Note that LM3-SNAP cannot simulate N₂O and NO emission rates. LM4.1-BNF underestimates dissolved organic N (DON) leaching rate (mean ∼0 vs. 0.6 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ leaching rate (mean ∼0 vs. 0.05 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in comparison to observations and reasonably estimates ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate (mean 0.1 vs. 0.1 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 6e) in comparison to observations. LM3-SNAP substantially overestimates ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate in comparison to observations (mean 8.3 vs. 0.1 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 6e).

Finally, LM4.1-BNF makes reasonable estimates for ecosystem C fluxes, particularly in comparison to LM3-SNAP. LM4.1-BNF effectively reproduces net ecosystem production (NEP) in both the growing and non-growing seasons at the hourly timescale, especially in comparison to LM3-SNAP, which overestimates NEP in the non-growing season (Fig. 7). LM4.1-BNF effectively reproduces gross primary production (GPP) and NPP, especially in comparison to LM3-SNAP, which substantially overestimates both GPP and NPP (Fig. 8). LM4.1-BNF overestimates GPP and NPP in comparison to observations (mean 15.4 vs. 13.1 and 8.1 vs. 7.5 $\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ respectively). LM4.1-BNF overestimates HR (mean 7.5 vs. 4.7 $\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and consequentially underestimates net ecosystem production (NEP) (mean 0.6 vs. 4.8 $\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in comparison to observations.

Figure 8 Simulated gross primary production (GPP), net primary production (NPP), heterotrophic respiration (HR), and net ecosystem production (NEP) by LM4.1-BNF and LM3-SNAP compared to the ForC database. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the data which are from mature forests. Error bars indicate 2 standard deviations.

4.2 Evaluation discussion

LM4.1-BNF effectively reproduces asymbiotic BNF and symbiotic BNF, i.e., its peak in early succession followed by its decline to zero in late succession (Fig. 3). In LM3-SNAP, asymbiotic BNF is not represented, and symbiotic BNF is sustained in late succession (Fig. 3). This occurs because there is no competitive exclusion of N-fixing plants by non-fixing plants in LM3-SNAP, which represents a single general plant C pool capable of BNF and cannot represent community dynamics (Fig. 4). In observed ecosystems, N-fixing plants are competitively excluded by non-fixing plants due to weak N limitation of plant growth in late succession (Menge et al., 2010; Sheffer et al., 2015), and symbiotic BNF is effectively zero (Boring and Swank, 1984).

This overestimation of symbiotic BNF in late succession by LM3-SNAP causes several problems. First, LM3-SNAP simulates lower total plant biomass C than LM4.1-BNF (mean 44 vs. 173 Mg C ha⁻¹ over the final 100 years; Fig. 5) because of the sustained high C cost of symbiotic BNF. Second, LM3-SNAP overestimates N losses in comparison to LM4.1-BNF, especially ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate (mean 8.3 vs. 0.1 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the final 100 years; Fig. 6). Accurately estimating N leaching rates is important given its downstream consequences such as eutrophication and acidification (Fowler et al., 2013; Tian and Niu, 2015).

5 LM4.1.-BNF performance relative to previous BNF representations

Here we describe the three analyses we conducted to identify LM4.1-BNF improvements to estimating major C and N pools and fluxes. First, we compare LM4.1-BNF to LM4.1-BNF with BNF represented as a function of NPP (LM4.1-BNF_NPP) and LM4.1-BNF with BNF represented as a function of ET (LM4.1-BNF_ET) to compare mechanistic and phenomenological representations of BNF. Second, we compare LM4.1-BNF simulations with both Robinia and Acer to LM4.1-BNF simulations with only Acer and LM4.1-BNF simulations with only Acer that can associate with N-fixing bacteria (hereafter, N-fixer Acer) to examine the role of competitive dynamics between N-fixing and non-fixing plants. Third, we compare LM4.1-BNF simulations with asymbiotic BNF to LM4.1-BNF simulations without asymbiotic BNF to examine the role of asymbiotic BNF.

5.1 Mechanistic and phenomenological representations of BNF

In our first analysis (Table 2), we compare mechanistic and phenomenological representations of BNF, and their implications for C and N cycling. LM4.1-BNF, LM4.1-BNF_NPP, and LM4.1-BNF_ET simulations estimate different total plant biomass C (Fig. 9a). LM4.1-BNF predicts the largest total plant biomass C (mean 170 Mg C ha⁻¹ over the final 100 years), followed by LM4.1-BNF_ET (mean 70 Mg C ha⁻¹ over the final 100 years) and LM4.1-BNF_NPP (mean 0 Mg C ha⁻¹ over the final 100 years). This is because, in LM4.1-BNF, BNF responds dynamically to strong N limitation of plant growth in early succession and BNF (mean 36 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the initial 100 years) supports total plant biomass C accumulation. Conversely, in LM4.1-BNF_NPP and LM4.1-BNF_ET, BNF does not respond dynamically to strong N limitation of plant growth in early succession, and BNF is not sufficient (mean 22 and ∼0 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the initial 100 years for LM4.1-BNF_ET and LM4.1-BNF_NPP respectively) to support total plant biomass C accumulation (Fig. 9b). As such, LM4.1-BNF effectively reproduces the temporal dynamics of symbiotic BNF rate, whereas LM4.1-BNF_ET and LM4.1-BNF_NPP predicted relatively constant symbiotic BNF rates. In observed ecosystems, strong N limitation of plant growth occurs in early succession. N-fixing trees are generally important pioneer species and can relieve strong N limitation of plant growth in early succession (Chapin III et al., 1994; Cierjacks et al., 2013; Menge et al., 2010). Consequently, symbiotic BNF is highest in early succession (Batterman et al., 2013; Boring and Swank, 1984; Menge and Hedin, 2009; Sullivan et al., 2014). Simulated soil C and N pools, soil N fluxes, soil N loss rates, and ecosystem C fluxes are relatively similar between simulations and are displayed in Figs. D5 and D6.

Figure 9 Simulated (a) total plant biomass C and (b) symbiotic BNF rate over time from LM4.1-BNF, LM4.1-BNF_NPP, and LM4.1-BNF_ET.

A similar result was found by Meyerholt et al. (2020), who compared five alternative representations of BNF within the O–CN model, including a BNF representation based on NPP, a BNF representation based on ET, and a BNF representation responding dynamically to N limitation of plant growth. As with our results, they found that the BNF representation responding dynamically to N limitation of plant growth predicted the largest total plant biomass C. However, their study did not compare these results to simulations that include competitive dynamics between N-fixing and non-fixing plants because O–CN represents a single general plant C pool capable of BNF and cannot represent community dynamics.

5.2 Competitive dynamics between N-fixing and non-fixing plants

In our second analysis (Table 2), we examine the role of competitive dynamics between N-fixing and non-fixing plants and its implication for C and N cycling. LM4.1-BNF simulations initialized with only Acer seedlings accumulate total plant biomass C slower than LM4.1-BNF simulations initialized with both Robinia and Acer seedlings (mean 48.6 vs. 80.9 Mg C ha⁻¹ over the initial 100 years respectively; Fig. 10a). In LM4.1-BNF simulations initialized with only Acer seedlings, stronger N limitation of plant growth in early succession due to the absence of a N-fixer vegetation type slows total plant biomass C accumulation. Nevertheless, total plant biomass C accumulates due to asymbiotic BNF and high N deposition at CHL (13.9 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), reaching a similar level to LM4.1-BNF simulations initialized with both Robinia and Acer seedlings after 100 years.

Figure 10 Simulated (a) total plant biomass C and (b) symbiotic BNF rate over time from LM4.1-BNF initialized with both Robinia and Acer, only Acer, and only N-fixer Acer, with and without asymbiotic BNF. aBNF indicates asymbiotic BNF.

LM4.1-BNF simulations initialized with only N-fixer Acer seedlings accumulate total plant biomass C similarly to LM4.1-BNF simulations initialized with both Robinia and Acer seedlings (mean 86.2 vs. 80.9 Mg C ha⁻¹ over the initial 100 years respectively; Fig. 10a). However, in LM4.1-BNF simulations initialized with only N-fixer Acer seedlings, a higher symbiotic BNF rate persists throughout succession in comparison to LM4.1-BNF simulations initialized with both Robinia and Acer seedlings (mean 23.6 vs. 0.2 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ over the final 100 years respectively; Fig. 10b). This occurs because there is no competitive exclusion of N-fixing plants by non-fixing plants due to weak N limitation of plant growth in late succession, which occurs in LM4.1-BNF simulations initialized with both Robinia and Acer seedlings. Simulated soil C and N pools, soil N fluxes, soil N loss rates, and ecosystem C fluxes are relatively similar between simulations and are displayed in Figs. D7 and D8.

Levy-Varon et al. (2019) conducted a similar study, in which a N-fixer vegetation type was included in the ED2 model. Similarly, they found that simulations without a N-fixer vegetation type accumulate total plant biomass C slower than simulations with a N-fixer vegetation type. However, ED2 differs from LM4.1-BNF in a multitude of processes. In particular, ED2 does not include representations of asymbiotic BNF, mycorrhizae, or rhizosphere priming. Furthermore, the representation of BNF in ED2 assumes instantaneous down-regulation of symbiotic BNF rate in comparison to the time lag in down-regulation of symbiotic BNF rate in LM4.1-BNF (due to the time between plant C allocation to symbiotic BNF, the growth of N-fixing bacteria, and symbiotic BNF) following observations (Bytnerowicz et al., 2021).

5.3 Asymbiotic BNF

In our third analysis (Table 2), we examine the role of asymbiotic BNF and its implications for C and N cycling. LM4.1-BNF simulations initialized with Acer without asymbiotic BNF accumulate total plant biomass C slower than LM4.1-BNF simulations initialized with Acer with asymbiotic BNF (30.4 vs. 48.6 Mg C ha⁻¹ over the initial 100 years; Fig. 10a). In LM4.1-BNF simulations initialized with Acer without asymbiotic BNF, stronger N limitation of plant growth in early succession due to the absence of both asymbiotic BNF and a N-fixer vegetation type (i.e., symbiotic BNF) substantially slows total plant biomass C accumulation. Nevertheless, total plant biomass C accumulates due to high N deposition at CHL (13.9 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), reaching a similar level to LM4.1-BNF simulations initialized with both Robinia and Acer seedlings and asymbiotic BNF after 300 years. Simulated soil C and N pools, soil N fluxes, soil N loss rates, and ecosystem C fluxes are relatively similar between simulations and are displayed in Figs. D7–8.

6 Discussion

6.1 Limitations

LM4.1-BNF captures the major pools and fluxes of C and N and their temporal dynamics. Importantly, LM4.1-BNF is novel in that it captures both the competitive dynamics between N-fixing and non-fixing plants and asymbiotic BNF. However, LM4.1-BNF has limitations.

LM4.1-BNF does not explicitly include asymbiotic BNF by bryophytes, lichens, and other organisms beyond soil microbes. This is regulated differently from asymbiotic BNF by soil microbes, specifically by light (Reed et al., 2011). In particular, in boreal forests and arctic tundra, asymbiotic BNF by bryophytes is a significant N flux (DeLuca et al., 2002). Additionally, herbaceous symbiotic BNF in the forest understory could be significant, but few studies have quantified its magnitude and controls (Cleveland et al., 1999).

The asymbiotic BNF temperature response is heavily biased towards high latitudes; the studies we used in its derivation had a mean latitude of 60^∘ (Chan, 1991; Chapin et al., 1991; Coxson and Kershaw, 1983; Liengen and Olsen, 1997; Roper, 1985). More studies on the asymbiotic BNF temperature response at lower latitudes are necessary.

The symbiotic BNF temperature response could acclimate to changing temperature (Bytnerowicz et al., 2021). The C cost of symbiotic BNF, which we assumed to be constant per unit N, could depend on temperature or other factors. These issues could influence the simulated response of symbiotic BNF and consequently total plant biomass C to increasing temperatures due to climate change. Thus, further empirical work on the effect of temperature on symbiotic BNF is necessary.

Finally, more observations of N cycling in general are necessary to validate N cycling representations in land models (Stocker et al., 2016; Thomas et al., 2015; Vicca et al., 2018). Global observations on N limitation of plant growth, soil N, N gas emission rates, N leaching rates, and, in particular, asymbiotic and symbiotic BNF rates are limited. Constraining these N pools and fluxes is critical to rigorously validating novel N cycling representations in land models.

6.2 Extensions

Robinia pseudoacacia is the most abundant N-fixing tree species in the coterminous United States (Staccone et al., 2020) and is also a common N-fixing tree across temperate forests; it is also found in Africa, Asia, Australia, Europe, and South America (Cierjacks et al., 2013). As such, Robinia pseudoacacia is representative of temperate N-fixing tree species.

The LM4.1-BNF representation of BNF, while implemented and validated in a temperate forest, can be expanded to other terrestrial ecosystems, such as tropical and boreal forests. This will require parameterization of representative N-fixing and non-fixing tree species but will not require re-structuring the model equations. Furthermore, the LM4.1-BNF representation of BNF could be incorporated into other land models.

Although N-fixing trees are generally important pioneer species and can relieve strong N limitation of plant growth in early succession (Chapin III et al., 1994), N-fixing trees can also be strong competitors. As such, in addition to having a facilitative effect on neighbouring plant growth (Hulvey et al., 2013), they can also have no effect on neighbouring plant growth (Lai et al., 2018; Xu et al., 2020) or an inhibitory effect on neighbouring plant growth (Chapin III et al., 2016; Taylor et al., 2017). This depends on abiotic and biotic factors (Staccone et al., 2021) and could be explored further with LM4.1-BNF.

7 Conclusions

Here we present LM4.1-BNF: an updated representation of BNF and other aspects of N cycling in LM4.1, which is the land component of the GFDL Earth System Model (Zhao et al., 2018a, b). LM4.1-BNF is the first land model to include a representation of the competitive dynamics between N-fixing and non-fixing plants, a mechanistic representation of asymbiotic BNF, and distinct asymbiotic and symbiotic BNF temperature responses derived from corresponding observations. Comparisons of simulations with observations show that LM4.1-BNF captures observed forest growth, successional dynamics, and major pools and fluxes of C and N and their temporal dynamics at population, community, and ecosystem scales. Furthermore, LM4.1-BNF represents these more accurately than previous representations of BNF in land models. By incorporating both the competitive dynamics between N-fixing and non-fixing plants, which is a key ecological mechanism that determines ecosystem-scale symbiotic BNF, and asymbiotic BNF, LM4.1-BNF yields accurate ecosystem-scale estimates of BNF and its temporal dynamics. Furthermore, the novel representation of soil NO and N₂O emissions in LM4.1-BNF enables the estimation of the magnitude of the terrestrial NO and N₂O source, which can be driven by BNF (Kou-Giesbrecht and Menge, 2021).

The representation of BNF in LM4.1-BNF is general and could be incorporated into other land models. Extending LM4.1-BNF to other biomes and incorporating LM4.1-BNF within the GFDL Earth System Model would allow a more accurate assessment of the response of BNF and the terrestrial C sink to elevated atmospheric CO₂ concentration, which intensifies N limitation of plant growth (Terrer et al., 2019; Zheng et al., 2020), and elevated N deposition, which relieves N limitation of plant growth (Reay et al., 2008; Schulte-Uebbing and de Vries, 2018; Zheng et al., 2020). In particular, such an endeavour could address whether BNF and N deposition will provide sufficient N to sustain CO₂ sequestration by terrestrial ecosystems under elevated atmospheric CO₂ concentration.

Appendix A: Model description

A1 N uptake by roots and symbionts

All vegetation types take up inorganic N via passive and active root uptake. Processes in Sect. A1 occur on the fast timescale (30 min).

A2 Passive root uptake of inorganic N

The rate of passive root uptake of inorganic N in soil layer k (N_passive(k); [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A1)}& {\mathrm{N}}_{\text{passive}}\left(k\right)=U_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(k\right)\left({\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)+{\mathrm{NH}}_{\mathrm{4}}\left(k\right)}{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\left(k\right)}}\right),\end{array}$$

where $U_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(k\right)$ is the water uptake flux in soil layer k [$\mathrm{kg}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$], NO₃(k) is the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool in soil layer k [kg N m⁻²], NH₄(k) is the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool in soil layer k [kg N m⁻²], and H₂O(k) is the soil water content of soil layer k [kg H₂O m⁻²]. Note that the solubility of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ differs, reflecting a higher cation exchange capacity than anion exchange capacity.

A2.1 Active root uptake of inorganic N

The rate of active root uptake of inorganic N in soil layer k (N_active(k); [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{N}}_{\text{active}}\left(k\right)=(r_{{\mathrm{NO}}_{\mathrm{3}}}{\displaystyle \frac{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)f_{\text{rhiz}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}}{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)f_{\text{rhiz}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)+k_{\mathrm{M},{\mathrm{NO}}_{\mathrm{3}}}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}}\\ \text{(A2)}& {\displaystyle }& {\displaystyle }+r_{{\mathrm{NH}}_{\mathrm{4}}}{\displaystyle \frac{\frac{{\mathrm{NH}}_{\mathrm{4}}\left(k\right)f_{\text{rhiz}}\left(k\right)}{\mathrm{\Delta }z}}{\frac{{\mathrm{NH}}_{\mathrm{4}}\left(k\right)f_{\text{rhiz}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}+k_{\mathrm{M},{\mathrm{NH}}_{\mathrm{4}}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}}){\mathrm{N}}_{\text{stress}},\end{array}$$

where $r_{{\mathrm{NO}}_{\mathrm{3}}}$ and $r_{{\mathrm{NH}}_{\mathrm{4}}}$ are rate constants, f_rhiz(k) is the rhizosphere volume fraction of soil layer k [m³ m⁻³], Δz(k) is the thickness of soil layer k, $k_{\mathrm{M},{\mathrm{NO}}_{\mathrm{3}}}$ and $k_{\mathrm{M},{\mathrm{NH}}_{\mathrm{4}}}$ are half-saturation constants, and n_indiv is the spatial density [indiv m⁻²]. f_rhiz(k) is calculated as

$$\begin{array}{}\text{(A3)}& f_{\text{rhiz}}\left(k\right)=\mathit{\pi }\left({\left(r_{\text{rhiz}}+r_{\text{root}}\right)}^{\mathrm{2}}-r_{\text{root}}^{\mathrm{2}}\right){\mathrm{C}}_{\mathrm{r}}\left(k\right)\text{SRL}{\displaystyle \frac{n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}},\end{array}$$

where r_rhiz is the radius of the rhizosphere around fine roots, r_root is the radius of fine roots, C_r(k) is the biomass C of fine roots in soil layer k [kg C indiv⁻¹], and SRL is the specific root length. N stress of the plant (N_stress; [unitless]) is given in Eq. (A34). This follows LM3-SNAP but is modified to increase with N stress.

A2.2 Inorganic N uptake by arbuscular mycorrhizae

The rates of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake by AM in soil layer k (${\mathrm{N}}_{\text{AM},{\mathrm{NO}}_{\mathrm{3}}}\left(k\right)$ and ${\mathrm{N}}_{\text{AM},{\mathrm{NH}}_{\mathrm{4}}}\left(k\right)$ respectively; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) are

$$\begin{array}{ll}{\displaystyle }{\mathrm{N}}_{\text{AM}}\left(k\right)=& {\displaystyle }{r}_{{\mathrm{NO}}_{\mathrm{3}},\text{AM}}{\displaystyle \frac{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}}{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}+k_{\text{AM},{\mathrm{NO}}_{\mathrm{3}}}}}\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{B_{\text{AM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}}{\frac{B_{\text{AM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}+k_{\text{AM}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}}\\ {\displaystyle }& {\displaystyle }+r_{{\mathrm{NH}}_{\mathrm{4}},\text{AM}}{\displaystyle \frac{\frac{{\mathrm{NH}}_{\mathrm{4}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}}{\frac{{\mathrm{NH}}_{\mathrm{4}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}+k_{\text{AM},{\mathrm{NH}}_{\mathrm{4}}}}}\\ \text{(A4)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{B_{\text{AM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}}{\frac{B_{\text{AM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}+k_{\text{AM}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}},\end{array}$$

where $r_{{\mathrm{NO}}_{\mathrm{3}},\text{AM}}$ and $r_{{\mathrm{NH}}_{\mathrm{4}},\text{AM}}$ are rate constants, $k_{\text{AM},{\mathrm{NO}}_{\mathrm{3}}}$ and $k_{\text{AM},{\mathrm{NH}}_{\mathrm{4}}}$ are half-saturation constants, B_AM(k) is the biomass C of AM in soil layer k [kg C indiv⁻¹], and k_AM is a half-saturation constant. This follows LM3-SNAP.

A2.3 Organic and inorganic N uptake by ectomycorrhizae

The rates of C and N decomposition by EM in soil layer k of organic matter type i, where i=labile plant-derived, labile microbe-derived, or recalcitrant ($D_{\mathrm{C},i,\text{EM}}\left(k\right)$ and $D_{\mathrm{N},i,\text{EM}}\left(k\right)$ respectively; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$] and [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]), are

$$\begin{array}{ll}{\displaystyle }{D}_{\mathrm{C},i,\text{EM}}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{EM},\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}}{\displaystyle \frac{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}{n_{\text{indiv}}}}\\ \text{(A5)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{B_{\text{EM}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{B_{\text{EM}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\text{M.EM}}}},{\displaystyle }{D}_{\mathrm{N},i,\text{EM}}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{EM},\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}}{\displaystyle \frac{{\mathrm{N}}_{\mathrm{U},i}\left(k\right)}{n_{\text{indiv}}}}\\ \text{(A6)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{B_{\text{EM}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{B_{\text{EM}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\text{M.EM}}}},\end{array}$$

where $V_{\text{EM},\text{max},\text{ref},i}$ is the maximum decomposition rate of organic matter type i, E_a,i is the activation energy of the decomposition of organic matter type i, R is the ideal gas constant, T_ref is reference temperature, T(k) is soil temperature of soil layer k [K], θ(k) is volumetric soil water content of soil layer k [m³ m⁻³], θ_sat is saturation volumetric soil water content [m³ m⁻³], f_θ,max is a factor normalizing the dependence on θ(k) to a maximum value of 1 [unitless], B_EM(k) is the biomass C of EM in soil layer k [kg C indiv⁻¹], C_U,i(k) is the soil C pool of type i in soil layer k [kg C m⁻²], N_U,i(k) is the soil N pool type i in soil layer k [kg C m⁻²], and k_M.EM is a half-saturation constant. This follows LM3-SNAP.

The rates of C and N uptake by EM in soil layer k (C_EM(k) and N_EM(k) respectively; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$] and [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) are

$$\begin{array}{ll}\text{(A7)}& {\displaystyle }{\mathrm{C}}_{\text{EM}}\left(k\right)=& {\displaystyle }\sum _i{\mathit{\epsilon }}_{\mathrm{C},i,\mathrm{EM}}{D}_{\mathrm{C},i,\mathrm{EM}}\left(k\right),{\displaystyle }{\mathrm{N}}_{\text{EM}}\left(k\right)=& {\displaystyle }\sum _i{\mathit{\epsilon }}_{\mathrm{N},i,\text{EM}}{D}_{\mathrm{N},i,\text{EM}}\left(k\right)+r_{{\mathrm{NO}}_{\mathrm{3}},\text{EM}}\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}}{\frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}+k_{\text{EM},{\mathrm{NO}}_{\mathrm{3}}}}}{\displaystyle \frac{\frac{B_{\text{EM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}}{\frac{B_{\text{EM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}+k_{\text{EM}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}}\\ {\displaystyle }& {\displaystyle }+r_{{\mathrm{NH}}_{\mathrm{4}},\text{EM}}{\displaystyle \frac{\frac{N{H}_{\mathrm{4}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}}{\frac{{\mathrm{NH}}_{\mathrm{4}}\left(k\right)}{\mathrm{\Delta }z\left(k\right)}+k_{\text{EM},{\mathrm{NH}}_{\mathrm{4}}}}}\\ \text{(A8)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{B_{\text{EM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}}{\frac{B_{\text{EM}}\left(k\right)n_{\text{indiv}}}{\mathrm{\Delta }z\left(k\right)}+k_{\text{EM}}}}{\displaystyle \frac{\mathrm{\Delta }z\left(k\right)}{n_{\text{indiv}}}},\end{array}$$

where ${\mathit{\epsilon }}_{\mathrm{C},i,\mathrm{EM}}$ is the C uptake efficiency of soil C type i by EM, ${\mathit{\epsilon }}_{\mathrm{N},i,\text{EM}}$ is the N uptake efficiency of soil N type i by EM, $r_{{\mathrm{NO}}_{\mathrm{3}},\text{EM}}$ and $r_{{\mathrm{NH}}_{\mathrm{4}},\text{EM}}$ are rate constants, $k_{\text{EM},{\mathrm{NO}}_{\mathrm{3}}}$ and $k_{\text{EM},{\mathrm{NH}}_{\mathrm{4}}}$ are half-saturation constants, and k_EM is a half-saturation constant. ${\sum }_i(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{C},i,\mathrm{EM}})D_{\mathrm{C},i,\mathrm{EM}}\left(k\right)$ is released as CO₂. ${\sum }_i(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{N},i,\text{EM}})D_{\mathrm{N},i,\text{EM}}\left(k\right)$ enters NH₄(k). This follows LM3-SNAP but is modified to additionally take up inorganic N.

A2.4 Symbiotic BNF by N-fixing bacteria

The symbiotic BNF rate by N-fixing bacteria (N_Nfix; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A9)}& {\mathrm{N}}_{\text{Nfix}}=r_{\text{Nfix}}{B}_{\text{Nfix}}f\left(T\right),\end{array}$$

where r_Nfix is a rate constant, B_Nfix is the biomass C of the nodule (includes both plant and N-fixing bacteria tissue) [kg C indiv⁻¹], and f_s(T) is the soil temperature dependence function. For Robinia,

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{f}_{\mathrm{s}}\left(T\right)=& max[\mathrm{0.0},\left({\displaystyle \frac{\mathrm{45.67}-\left(T-\mathrm{273.15}\right)}{\mathrm{45.67}-\mathrm{31.89}}}\right)\\ & {\left({\displaystyle \frac{\left(T-\mathrm{273.15}\right)-\mathrm{1.43}}{\mathrm{31.89}-\mathrm{1.43}}}\right)}^{\frac{\mathrm{31.89}-\mathrm{1.43}}{\mathrm{45.67}-\mathrm{31.89}}}],\end{array}\end{array}$$

where T is the average soil temperature across soil layers [K]. This reaches its maximum at 31.9 ^∘C (Fig. 2). This is derived from Bytnerowicz et al. (2021).

A3 Asymbiotic BNF

Soil microbes are represented as a single C pool that conducts decomposition, nitrification, denitrification, and asymbiotic BNF. The rates of C and N decomposition by soil microbes in soil layer k of organic matter type i, where i=labile plant-derived, labile microbe-derived, or recalcitrant, (D_C,i(k) and D_N,i(k) respectively; [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] and [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]), are

$$\begin{array}{ll}{\displaystyle }{D}_{\mathrm{C},i}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}}{\mathrm{C}}_{\mathrm{U},i}\left(k\right)\\ \text{(A11)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\mathrm{M}}}},{\displaystyle }{D}_{\mathrm{N},i}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}}{\mathrm{N}}_{\mathrm{U},i}\left(k\right)\\ \text{(A12)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\mathrm{M}}}},\end{array}$$

where $V_{\text{max},\text{ref},i}$ is the maximum decomposition rate of organic matter type i, C_M(k) is the biomass C of soil microbes in soil layer k [kg C m⁻²], and k_M is the half-saturation constant. This follows LM3-SNAP.

The potential rates of C and N decomposition during denitrification by soil microbes in soil layer k of organic matter type i ($D_{\mathrm{C},i,\text{denit},\text{pot}}\left(k\right)$ and $D_{\mathrm{N},i,\text{denit},\text{pot}}\left(k\right)$ respectively; [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] and [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) are

$$\begin{array}{ll}{\displaystyle }{D}_{\mathrm{C},i,\text{denit},\text{pot}}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{denit},\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{5.5}}}{f_{\mathit{\theta },\text{max}}}}{\mathrm{C}}_{\mathrm{U},i}\left(k\right)\\ \text{(A13)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\mathrm{M},\text{denit}}}},{\displaystyle }{D}_{\mathrm{N},i,\text{denit},\text{pot}}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{denit},\text{max},\text{ref},i}}{\exp \left(\frac{-E_{\mathrm{a},i}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},i}}{RT\left(k\right)}}\right)\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{5.5}}}{f_{\mathit{\theta },\text{max}}}}{\mathrm{N}}_{\mathrm{U},i}\left(k\right)\\ \text{(A14)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}}{\frac{C_{\mathrm{M}}\left(k\right)}{{\mathrm{C}}_{\mathrm{U},i}\left(k\right)}+k_{\mathrm{M},\text{denit}}}},\end{array}$$

where $V_{\text{denit},\text{max},\text{ref},i}$ is the maximum decomposition rate of organic matter type i during denitrification and k_M,denit is the half-saturation constant. This follows LM3-SNAP.

The rates of C and N decomposition during denitrification by soil microbes in soil layer k of organic matter type i ($D_{\mathrm{C},i,\text{denit}}\left(k\right)$ and $D_{\mathrm{N},i,\text{denit}}\left(k\right)$ respectively; [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] and [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) are

$$\begin{array}{}\text{(A15)}& {\displaystyle }{D}_{\mathrm{C},i,\text{denit}}\left(k\right)={\displaystyle \frac{D_{\mathrm{C},i,\text{denit},\text{pot}}\left(k\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)+k_{\text{denit}}{f}_{\text{denit}}{\sum }_i{D}_{\mathrm{C},i,\text{denit},\text{pot}}\left(k\right)}},\text{(A16)}& {\displaystyle }{D}_{\mathrm{N},i,\text{denit}}\left(k\right)={\displaystyle \frac{D_{\mathrm{N},i,\text{denit},\text{pot}}\left(k\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)+k_{\text{denit}}{f}_{\text{denit}}{\sum }_i{D}_{\mathrm{C},i,\text{denit},\text{pot}}\left(k\right)}},\end{array}$$

where k_denit is the half-saturation constant, and f_denit is the stoichiometric ratio of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ demand for C decomposition. $D_{\mathrm{C},i,\text{denit}}\left(k\right)$ and $D_{\mathrm{N},i,\text{denit}}\left(k\right)$ include ${\mathrm{NO}}_{\mathrm{3}}^{-}$ limitation of denitrification. This follows LM3-SNAP.

The rate of change of biomass C of soil microbes in soil layer k ($\left(\mathrm{d}{\mathrm{C}}_{\mathrm{M}}\right(k\left)\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A17)}& {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{M}}\left(k\right)}{\mathrm{d}t}}=G_{\mathrm{M},\mathrm{C}}\left(k\right)-{\displaystyle \frac{{\mathrm{C}}_{\mathrm{M}}\left(k\right)}{{\mathit{\tau }}_{\mathrm{M}}}},\end{array}$$

where G_M,C(k) is the C growth rate of soil microbes in soil layer k (Eq. A23) and τ_M is the combined maintenance respiration and turnover time of soil microbes. This follows LM3-SNAP.

The rate of change of biomass N of soil microbes in soil layer k ($\left(\mathrm{d}{\mathrm{N}}_{\mathrm{M}}\right(k\left)\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A18)}& {\displaystyle \frac{\mathrm{d}{\mathrm{N}}_{\mathrm{M}}\left(k\right)}{\mathrm{d}t}}=G_{\mathrm{M},\mathrm{N}}\left(k\right)-{\displaystyle \frac{\frac{{\mathrm{C}}_{\mathrm{M}}\left(k\right)}{{\mathit{\tau }}_{\mathrm{M}}}{\mathit{\epsilon }}_{t,\mathrm{M}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{M}}}},\end{array}$$

where G_M,N(k) is the N growth rate of soil microbes in soil layer k (Eq. A24), ε_t,M is the fraction of maintenance respiration in combined maintenance respiration and turnover, and C:N_M is the C:N ratio of soil microbes. This follows LM3-SNAP.

The maintenance respiration rate of soil microbes in soil layer k (R_maint(k); [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A19)}& R_{\text{maint}}\left(k\right)={\displaystyle \frac{{\mathrm{C}}_{\mathrm{M}}\left(k\right)}{{\mathit{\tau }}_{\mathrm{M}}}}\left(\mathrm{1}-{\mathit{\epsilon }}_{t,\mathrm{M}}\right).\end{array}$$

This is released as CO₂. This follows LM3-SNAP.

The asymbiotic BNF rate of soil microbes in soil layer k (N_{Nfix asymb}(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A20)}& {\mathrm{N}}_{\text{Nfix asymb}}\left(k\right)=r_{\text{Nfix asymb}}{\mathrm{C}}_{\mathrm{M}}\left(k\right)f\left(T\left(k\right)\right),\end{array}$$

where r_{Nfix asymb} is a rate constant and f_a(T(k)) is the soil temperature dependence function.

$$\begin{array}{}\text{(A21)}& f_a\left(T\left(k\right)\right)=e^{-\mathrm{2.6}+\mathrm{0.21}\left(T\left(k\right)-\mathrm{273.15}\right)\left(\mathrm{1}-\frac{\mathrm{0.5}\left(T\left(k\right)-\mathrm{273.15}\right)}{\mathrm{24.4}}\right)},\end{array}$$

which reaches its maximum at 24.4 ^∘C (Fig. 2). This is derived from the observations compiled by Houlton et al. (2008) with the study of symbiotic BNF removed (Schomberg and Weaver, 1992) and is normalized to a maximum of 1.

The N surplus or deficit of soil microbes in soil layer k (ϕ_N(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_{\mathrm{N}}\left(k\right)=& {\displaystyle }\sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i}\left(k\right)+\sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i,\text{denit}}\left(k\right)\\ {\displaystyle }& {\displaystyle }+{\mathrm{N}}_{\text{Nfix}\text{asymb}}\left(k\right)\\ \text{(A22)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\sum }_i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i}\left(k\right)+{\sum }_i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i,\text{denit}}\left(k\right)-R_{\text{maint}}\left(k\right)}{{\mathrm{C}:\mathrm{N}}_{\mathrm{M}}}},\end{array}$$

where ε_N,i is the N uptake efficiency of soil N type i by soil microbes, and ε_C,i is the C uptake efficiency of soil C type i by soil microbes. ϕ_N(k)>0 indicates net N mineralization (N surplus), and ϕ_N(k)<0 indicates net N immobilization (N deficit) by soil microbes in soil layer k. ${\sum }_i(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{C},i})\left(D_{\mathrm{C},i}\right(k)+D_{\mathrm{C},i,\text{denit}}(k\left)\right)$ is released as CO₂. ${\sum }_i(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{N}})\left(D_{\mathrm{N},i}\right(k)+D_{\mathrm{N},i,\text{denit}}(k\left)\right)$ enters NH₄(k). This follows LM3-SNAP but was modified to include asymbiotic BNF.

G_M,C(k) and G_M,N(k) depend on whether growth of soil microbes is C-limited (${\mathit{\varphi }}_{\mathrm{N}}\left(k\right)\ge -{\text{Imm}}_{\text{max}}\left(k\right))$ or N-limited (${\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right))$ and are calculated as follows.

$$\begin{array}{}\text{(A23)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{M},\mathrm{C}}\left(k\right)=\left\{\begin{array}{l}\sum _i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i}\left(k\right)\\ +\sum _i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i,\text{denit}}\left(k\right),& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)\ge -{\text{Imm}}_{\text{max}}\left(k\right)\\ {\mathrm{C}:\mathrm{N}}_{\mathrm{M}}\left[\sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i}\right(k)\\ +\sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i,\text{denit}}\left(k\right)\\ +{\mathrm{N}}_{\text{Nfix},\text{asymb}}\left(k\right)\\ +{\text{Imm}}_{\text{max}}\left(k\right)]\\ +R_{\text{maint}}\left(k\right),& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right)\end{array}\right.\text{(A24)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{M},\mathrm{N}}\left(k\right)=\left\{\begin{array}{l}\left(\sum _i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i}\right(k)\\ +\sum _i{\mathit{\epsilon }}_{\mathrm{C},i}{D}_{\mathrm{C},i,\text{denit}}\left(k\right)\\ -R_{\text{maint}}\left(k\right))\\ /\left({\mathrm{C}:\mathrm{N}}_{\mathrm{M}}\right),& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)\ge -{\text{Imm}}_{\text{max}}\left(k\right)\\ \sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i}\left(k\right)\\ +\sum _i{\mathit{\epsilon }}_{\mathrm{N},i}{D}_{\mathrm{N},i,\text{denit}}\left(k\right)\\ +{\mathrm{N}}_{\text{Nfix},\text{asymb}}\left(k\right)\\ +{\text{Imm}}_{\text{max}}\left(k\right),& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right)\end{array}\right.\end{array}$$

This follows LM3-SNAP but was modified to include asymbiotic BNF.

The maximum N immobilization rate of soil microbes in soil layer k (Imm_max(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }{\text{Imm}}_{\text{max}}\left(k\right)=& {\displaystyle }\left(V_{\text{max},\text{ref},{\mathrm{NH}}_{\mathrm{4}}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{RT\left(k\right)}}\right){\mathrm{NH}}_{\mathrm{4}}\right(k)\\ {\displaystyle }& {\displaystyle }+V_{\text{max},\text{ref},{\mathrm{NO}}_{\mathrm{3}}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{RT\left(k\right)}}\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right))\\ \text{(A25)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}},\end{array}$$

where $V_{\text{max},\text{ref},{\mathrm{NH}}_{\mathrm{4}}}$ is the maximum ${\mathrm{NH}}_{\mathrm{4}}^{+}$ immobilization rate, $E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}$ is the activation energy of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ immobilization, $V_{\text{max},\text{ref},{\mathrm{NO}}_{\mathrm{3}}}$ is the maximum ${\mathrm{NO}}_{\mathrm{3}}^{-}$ immobilization rate, and $E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}$ is the activation energy of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ immobilization. This follows LM3-SNAP.

When growth of soil microbes in soil layer k is N-limited (${\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right))$, there is overflow respiration of excess C. The overflow respiration of excess C in soil layer k (R_overflow(k); [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A26)}& R_{\text{overflow}}\left(k\right)=\left\{\begin{array}{ll}\mathrm{0},& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)\ge -{\text{Imm}}_{\text{max}}\left(k\right)\\ -\left({\mathit{\varphi }}_{\mathrm{N}}\right(k)\\ +{\text{Imm}}_{\text{max}}\left(k\right))\\ \cdot \left({\mathrm{C}:\mathrm{N}}_{\mathrm{M}}\right),& {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right)\end{array}\right..\end{array}$$

This is released as CO₂.

If ϕ_N(k)>0, the net N mineralization flux in soil layer k is ϕ_N(k). If $-{\text{Imm}}_{\text{max}}\left(k\right)\le {\mathit{\varphi }}_{\mathrm{N}}\left(k\right)\le \mathrm{0}$, the net N immobilization flux in soil layer k is −ϕ_N(k). If ${\mathit{\varphi }}_{\mathrm{N}}\left(k\right)<-{\text{Imm}}_{\text{max}}\left(k\right)$, the net N immobilization flux in soil layer k is −Imm_max(k). The fraction of the net N immobilization flux in soil layer k that is ${\mathrm{NH}}_{\mathrm{4}}^{+}$ immobilization is

$${\displaystyle \frac{\frac{V_{\text{max},\text{ref},{\mathrm{NH}}_{\mathrm{4}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{RT\left(k\right)}\right){\mathrm{NH}}_{\mathrm{4}}\left(k\right)}{\frac{V_{\text{max},\text{ref},{\mathrm{NH}}_{\mathrm{4}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{RT\left(k\right)}\right){\mathrm{NH}}_{\mathrm{4}}\left(k\right)+\frac{V_{\text{max},\text{ref},{\mathrm{NO}}_{\mathrm{3}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{RT\left(k\right)}\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right)}}.$$

The fraction of the net N immobilization flux in soil layer k that is ${\mathrm{NO}}_{\mathrm{3}}^{-}$ immobilization is

$${\displaystyle \frac{\frac{V_{\text{max},\text{ref},{\mathrm{NO}}_{\mathrm{3}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{RT\left(k\right)}\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\frac{V_{\text{max},\text{ref},{\mathrm{NH}}_{\mathrm{4}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NH}}_{\mathrm{4}}}}{RT\left(k\right)}\right){\mathrm{NH}}_{\mathrm{4}}\left(k\right)+\frac{V_{\text{max},\text{ref},{\mathrm{NO}}_{\mathrm{3}}}}{\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{R{T}_{\text{ref}}}\right)}\exp \left(\frac{-E_{\mathrm{a},{\mathrm{NO}}_{\mathrm{3}}}}{RT\left(k\right)}\right){\mathrm{NO}}_{\mathrm{3}}\left(k\right)}}.$$

This follows LM3-SNAP.

Processes in Sect. A2 occur on the fast timescale (30 min).

A4 Plant growth and N limitation

The non-structural C pool (NSC; [kg C indiv⁻¹]) gains C from photosynthesis. NSC loses C to respiration and C allocation to growth, symbionts, and root C exudation. The rate of change of NSC ($\left(\text{dNSC}\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle \frac{\text{dNSC}}{\mathrm{d}t}}=& {\displaystyle }P-R-\left(G_{\mathrm{C},\mathrm{l}}+G_{\mathrm{C},\mathrm{r}}+G_{\mathrm{C},\text{sw}}+G_{\mathrm{C},\text{seed}}\right)\\ \text{(A27)}& {\displaystyle }& {\displaystyle }-{\mathrm{C}}_{\text{alloc}}-L_{\mathrm{C},\text{exudate}},\end{array}$$

where P is the photosynthesis rate [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], R is the respiration rate (maintenance and growth) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,l is the growth rate of the leaf C pool (C_l; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,r is the growth rate of the fine-root C pool (C_r; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,sw is the growth rate of the sapwood C pool (C_sw; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], G_C,seed is the growth rate of the seed C pool (C_seed; [kg C indiv⁻¹]) [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], C_alloc is the rate of C allocation to symbionts (Eqs. A43–A46), and L_C,exudate is the rate of root C exudation (Eq. A38). Note that sapwood is converted to heartwood following Martinez Cano et al. (2020).

The non-structural N pool (NSN; [kg N indiv⁻¹]) gains N from N uptake via roots and symbionts. NSN loses N to N allocation to growth, symbionts, and root N exudation. The rate of change of NSN ($\left(\text{dNSN}\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle \frac{\text{dNSN}}{\mathrm{d}t}}=& {\displaystyle }U-\left({\displaystyle \frac{G_{\mathrm{C},\mathrm{l}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{l}}}}+{\displaystyle \frac{G_{\mathrm{C},\mathrm{r}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{r}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{sw}}}{{\mathrm{C}:\mathrm{N}}_{\text{sw}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{seed}}}{{\mathrm{C}:\mathrm{N}}_{\text{seed}}}}\right)\\ \text{(A28)}& {\displaystyle }& {\displaystyle }-{\mathrm{N}}_{\text{alloc}}-L_{\mathrm{N},\text{exudate}},\end{array}$$

where U is the N uptake rate via roots and symbionts (Eq. A55) [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$], C:N_l is the fixed C:N ratio of leaves, C:N_r is the fixed C:N ratio of fine roots, C:N_sw is the fixed C:N ratio of sapwood, C:N_seed is the fixed C:N ratio of seeds, N_alloc is the rate of N allocation to symbionts (Eq. A47), and L_N,exudate is the rate of root N exudation (Eq. A39).

Non-N-limited growth is calculated according to Weng et al. (2015). The total allocation of NSC to growth is determined by the target NSC (NSC_target; [kg C indiv⁻¹]) and minimizes the deviation between NSC and NSC_target. NSC_target is a multiple of the target C_l (C_l,target; [kg C indiv⁻¹]), which reflects the ability of a plant to refoliate after defoliation (Hoch et al., 2003; Richardson et al., 2013), and is calculated as

$$\begin{array}{}\text{(A29)}& {\text{NSC}}_{\text{target}}=q{\mathrm{C}}_{\mathrm{l},\text{target}},\end{array}$$

where q is a proportionality constant. The allocation of NSC to the growth of each tissue depends on the total allocation of NSC to growth and the target C pool of each tissue and minimizes the deviation between the C pool of each tissue and the target C pool of each tissue. The target C pool of each tissue is dynamic and is determined by allometry (Eqs. A40–A42).

In LM4.1-BNF, G_C,l, G_C,r, G_C,sw, and G_C,seed are adjusted to include N limitation and are calculated as

$$\begin{array}{}\text{(A30)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\mathrm{l}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\mathrm{l}},\text{(A31)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\mathrm{r}}={\mathrm{\Delta }}_{\mathrm{r}},\text{(A32)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\text{sw}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\text{sw}},\text{(A33)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{C},\text{seed}}=(\mathrm{1}-{\mathrm{N}}_{\text{stress}}){\mathrm{\Delta }}_{\text{seed}},\end{array}$$

where N_stress is N stress [unitless] and Δ_l, Δ_r, Δ_sw, and Δ_seed are the non-N-limited growth rates of C_l, C_r, C_sw, and C_seed respectively [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$] following Weng et al. (2015). Because plants increase C allocation to fine roots relative to other tissues when N-limited (Poorter et al., 2012), G_C,r is not adjusted to include N limitation.

In LM4.1-BNF, N_stress is the relative difference between NSN and NSN_target and is calculated as

$$\begin{array}{}\text{(A34)}& {\mathrm{N}}_{\text{stress}}=max\left[\mathrm{0},{\displaystyle \frac{{\text{NSN}}_{\text{target}}-\text{NSN}}{{\text{NSN}}_{\text{target}}}}\right],\end{array}$$

where NSN_target is the target NSN [kg N indiv⁻¹]. N_stress is smoothed with a low-pass filter over 30 d to reflect the persisting influence of N stress (Mooney et al., 1991). NSN_target is calculated as

$$\begin{array}{}\text{(A35)}& {\text{NSN}}_{\text{target}}={\displaystyle \frac{{\text{NSC}}_{\text{target}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{l}}}}.\end{array}$$

This is similar to LM3-SNAP, which compared the target leaf and root N pools to NSN, but is modified to reflect the treatment of NSC_target in LM4.1 by including the target sapwood and seed N pools.

N demand for plant growth (N_demand; [kg N indiv⁻¹]) is

$$\begin{array}{}\text{(A36)}& {\mathrm{N}}_{\text{demand}}={\displaystyle \frac{G_{\mathrm{C},\mathrm{l}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{l}}}}+{\displaystyle \frac{G_{\mathrm{C},\mathrm{r}}}{{\mathrm{C}:\mathrm{N}}_{\mathrm{r}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{sw}}}{{\mathrm{C}:\mathrm{N}}_{\text{sw}}}}+{\displaystyle \frac{G_{\mathrm{C},\text{seed}}}{{\mathrm{C}:\mathrm{N}}_{\text{seed}}}}.\end{array}$$

If N_demand>0.5NSN, $G_{\mathrm{C},\mathrm{l}},G_{\mathrm{C},\mathrm{r}},G_{\mathrm{C},\text{sw}}$, and G_C,seed are reduced to prevent N_demand from depleting NSN

$$\begin{array}{}\text{(A37)}& G_{\mathrm{C},t}=\left\{\begin{array}{ll}{G}_{\mathrm{C},t},& {\mathrm{N}}_{\text{demand}}<\mathrm{0.5}\text{NSN}\\ G_{\mathrm{C},t}\frac{\text{NSN}}{{\mathrm{N}}_{\text{demand}}},& {\mathrm{N}}_{\text{demand}}>\mathrm{0.5}\text{NSN},\end{array}\right.\end{array}$$

where t=leaf, root, sapwood, or seed, and 0.5 was set to maintain a baseline NSN.

Plant turnover decreases C_l, C_r, and C_sw and from N_l, N_r, and N_sw at a constant tissue-specific rate and enters C_{U,labile plant-derived} or C_{U,recalcitrant} and N_{U,labile plant-derived} or N_{U,recalcitrant} respectively. A fraction of the turnover of C_l and N_l is retranslocated into NSC and NSN respectively.

Under N limitation, plants increase root C exudation to stimulate N mineralization in the rhizosphere (rhizosphere priming; Cheng et al., 2014; Finzi et al., 2015). L_C,exudate increases with N_stress and is calculated as

$$\begin{array}{}\text{(A38)}& L_{\mathrm{C},\text{exudate}}=r_{\text{leakage},\mathrm{C}}\text{NSC}{\mathrm{N}}_{\text{stress}},\end{array}$$

where r_leakage,C is a rate constant. L_C,exudate enters the rhizosphere C_{U,labile plant-derived}.

Under N limitation, plants decrease root N exudation (Canarini et al., 2019). L_N,exudate decreases with N_stress and is calculated as

$$\begin{array}{}\text{(A39)}& L_{\mathrm{N},\text{exudate}}=r_{\text{leakage},\mathrm{N}}\text{NSN}\left(\mathrm{1}-{\mathrm{N}}_{\text{stress}}\right),\end{array}$$

where r_leakage,N is a rate constant. L_N,exudate enters the rhizosphere N_{U,labile plant-derived}.

Plant growth occurs during the growing season. Plant maintenance respiration and plant turnover occur throughout the year. After the transition from the growing season to the non-growing season, turnover of C_l and N_l occur at an elevated rate until C_l=0 and N_l=0. Plant C allocation to symbionts occurs during the growing season. Root exudation occurs during the growing season. Symbiont growth, maintenance respiration, and turnover occur throughout the year. N uptake by roots and symbionts occurs throughout the year.

Photosynthesis and respiration occur on the fast timescale (30 min). Plant turnover occurs on the fast timescale (30 min). Plant growth, plant C allocation to symbionts, and root exudation occur on the daily timescale.

A5 Plant allometry

Plant allometry follows Martinez Cano et al., 2020. Crown area (CA; [m²]) is a function of diameter at breast height (power function):

$$\begin{array}{}\text{(A40)}& \text{CA}={\mathit{\alpha }}_{\text{CA}}{D}^{{\mathit{\theta }}_{\text{CA}}},\end{array}$$

where D is diameter at breast height [m], and α_CA and θ_CA are allometry parameters.

Height (H; [m]) is a function of diameter at breast height (generalized Michaelis–Menten equation)

$$\begin{array}{}\text{(A41)}& H={\displaystyle \frac{{\mathit{\alpha }}_{\text{HT}}{D}^{{\mathit{\theta }}_{\text{HT}}}}{{\mathit{\gamma }}_{\text{HT}}+D^{{\mathit{\theta }}_{\text{HT}}}}},\end{array}$$

where α_HT, θ_HT, and γ_HT are allometry parameters.

Wood mass (W; [kg C]) is a function of diameter at breast height and height:

$$\begin{array}{}\text{(A42)}& W={\mathit{\alpha }}_{\text{BM}}{\mathit{\rho }}_{\text{wood}}{D}^{\mathrm{2}}H,\end{array}$$

where α_BM is an allometry parameter and ρ_wood is wood density.

A6 Plant C allocation to symbionts (AM, EM, and N-fixing bacteria)

The rate of C allocation to AM (C_alloc,AM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A43)}& {\mathrm{C}}_{\text{alloc},\text{AM}}=f_{\text{alloc},\text{AM}}\text{NSC},\end{array}$$

where f_alloc,AM is the fraction of NSC allocated to AM per unit time. C_alloc,AM is not related to N_stress because, although AM increases N uptake, AM is maintained by the plant primarily for phosphorus uptake (Smith and Smith, 2011).

The rate of C allocation to EM (C_alloc,EM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A44)}& {\mathrm{C}}_{\text{alloc},\text{EM}}=f_{\text{alloc},\text{EM}}\text{NSC}{\mathrm{N}}_{\text{stress}},\end{array}$$

where f_alloc,EM is the maximum fraction of NSC allocated to EM per unit time. C_alloc,EM is a function of N_stress because biomass C of EM increases with N limitation (Phillips et al., 2013).

Plants that associate with N-fixing bacteria can regulate symbiotic BNF to different extents, termed their BNF strategy (Menge et al., 2015). For plants with a perfectly facultative BNF strategy, symbiotic BNF increases with N limitation. For plants with an incomplete BNF strategy, symbiotic BNF increases with N limitation but is maintained at a minimum. For plants with an obligate BNF strategy, symbiotic BNF is constant. For plants with either a facultative or an incomplete BNF strategy, the rate of C allocation by the plant to N-fixing bacteria (C_alloc,Nfix; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A45)}& {\mathrm{C}}_{\text{alloc},\text{Nfix}}=max\left[f_{\text{alloc},\text{Nfix}}\text{NSC}{\mathrm{N}}_{\text{stress}},f_{\text{alloc},\text{Nfix},\text{min}}\text{NSC}\right],\end{array}$$

where f_alloc,Nfix is the fraction of NSC allocated to N-fixing bacteria per unit time and $f_{\text{alloc},\text{Nfix},\text{min}}$ is the minimum fraction of NSC allocated to N-fixing bacteria per unit time. For plants with a perfectly facultative BNF strategy $f_{\text{alloc},\text{Nfix},\text{min}}=\mathrm{0}$, and for plants with an incomplete down-regulator BNF strategy $f_{\text{alloc},\text{Nfix},\text{min}}>\mathrm{0}$. We model Robinia with an incomplete down-regulator BNF strategy.

For plants with an obligate BNF strategy, C_alloc,Nfix is

$$\begin{array}{}\text{(A46)}& {\mathrm{C}}_{\text{alloc},\text{Nfix}}=f_{\text{alloc},\text{Nfix}}\text{NSC}.\end{array}$$

Additionally, plants allocate a small quantity of N to symbionts such that symbiont growth can be initiated. The rate of N allocation by the plant to symbionts (N_alloc,j; $j=\text{AM},$ EM, N-fixing bacteria; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A47)}& {\mathrm{N}}_{\text{alloc},j}={\displaystyle \frac{{\mathrm{C}}_{\text{alloc},j}}{{\mathrm{C}:\mathrm{N}}_{\text{alloc}}}},\end{array}$$

where C:N_alloc is the C:N ratio of C and N allocated to symbionts by the plant.

Processes in Sect. A5 occur on the daily timescale.

A7 Growth and turnover of symbionts

Plant C allocation to symbionts is transferred to an intermediate C pool (C_int,j; j=AM, EM, or N-fixing bacteria; [kg C indiv⁻¹]). The rate of change of C_int,j ($\left(\mathrm{d}{\mathrm{C}}_{\text{int},j}\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A48)}& {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\text{int},j}}{\mathrm{d}t}}={\mathrm{C}}_{\text{alloc},j}-{\displaystyle \frac{G_j}{{\mathit{\epsilon }}_{\text{symb}}}},\end{array}$$

where ε_symb is the proportion of C uptake by a symbiont from C_int,j that is assimilated. The C growth rate of a symbiont with strategy j (G_j; j=AM or N-fixing bacteria; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A49)}& G_j={\mathit{\epsilon }}_{\text{symb}}{r}_{\text{growth}}{\mathrm{C}}_{\text{int},j},\end{array}$$

where r_growth is the growth rate of a symbiont.

The growth rate of EM (G_EM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A50)}& G_{\text{EM}}={\mathit{\epsilon }}_{\text{symb}}(r_{\text{growth}}{\mathrm{C}}_{\text{int},\text{EM}}+{\mathrm{C}}_{\text{EM}}).\end{array}$$

$\left(G_j\right)/\left({\mathit{\epsilon }}_{\text{symb}}\right)(\mathrm{1}-{\mathit{\epsilon }}_{\text{symb}})$ (j=AM, EM, or N-fixing bacteria) is released as CO₂ (growth respiration).

The rate of change of biomass C of a symbiont with strategy j ($\left(\mathrm{d}{B}_j\right)/\left(\mathrm{d}t\right)$; j=AM or EM; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A51)}& {\displaystyle \frac{\mathrm{d}{B}_j}{\mathrm{d}t}}=G_j-{\mathit{\xi }}_j{B}_j-{\displaystyle \frac{B_j}{{\mathit{\tau }}_j}},\end{array}$$

where ξ_j is rate of maintenance respiration of a symbiont with strategy j, and τ_j is the turnover time of a symbiont with strategy j. ξ_jB_j is released as CO₂ (maintenance respiration). $\left(B_j\right)/\left({\mathit{\tau }}_j\right)$ enters C_{U,labile microbe-derived} and $\left(B_j\right)/\left({\mathrm{C}:\mathrm{N}}_j{\mathit{\tau }}_j\right)$ enters N_{U,labile microbe-derived}, where C:N_j is the C:N ratio of a symbiont with strategy j.

The rate of change of biomass C of N-fixing bacteria ($\left(\mathrm{d}{B}_{\text{Nfix}}\right)/\left(\mathrm{d}t\right)$; [$\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A52)}& {\displaystyle \frac{\mathrm{d}{B}_{\text{Nfix}}}{\mathrm{d}t}}=G_{\text{Nfix}}-{\mathit{\xi }}_{\text{Nfix}}{B}_{\text{Nfix}}-{\text{cost}}_{\text{Nfix}}{\mathrm{N}}_{\text{Nfix}}-{\displaystyle \frac{B_{\text{Nfix}}}{{\mathit{\tau }}_{\text{Nfix}}}},\end{array}$$

where ξ_Nfix is rate of maintenance respiration of N-fixing bacteria, cost_Nfix is the C cost of symbiotic BNF per unit N, and τ_Nfix is the turnover time of N-fixing bacteria. ξ_NfixB_Nfix is released as CO₂ (maintenance respiration), and cost_NfixN_Nfix is released as CO₂ (respiration associated with symbiotic BNF). $\left(B_{\text{Nfix}}\right)/\left({\mathit{\tau }}_{\text{Nfix}}\right)$ enters C_{U,labile microbe-derived}, and $\left(B_{\text{Nfix}}\right)/\left({\mathrm{C}:\mathrm{N}}_{\text{Nfix}}{\mathit{\tau }}_{\text{Nfix}}\right)$ enters N_{U,labile microbe-derived}, where C:N_Nfix is the C:N ratio of N-fixing bacteria.

N acquired by symbionts is transferred to an intermediate N pool (N_int,j; [kg N indiv⁻¹]). The rate of change of N_int,j ($\left(\mathrm{d}{\mathrm{N}}_{\text{int},j}\right)/\left(\mathrm{d}t\right)$; j=AM or EM; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A53)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{N}}_{\text{int},j}}{\mathrm{d}t}}& =\sum _k{\mathrm{N}}_j+{\displaystyle \frac{{\mathrm{C}}_{\text{alloc},j}}{{\mathrm{C}:\mathrm{N}}_{\text{alloc}}}}-{\displaystyle \frac{\mathrm{1}}{{\mathrm{C}:\mathrm{N}}_j}}\left(G_j-{\mathit{\xi }}_j{B}_j\right)\\ & -r_{\text{up},\text{veg}}{\mathrm{N}}_{\text{int},j},\end{array}\end{array}$$

where r_up,veg is the rate of plant N uptake from N_int,j.

The rate of change of N_int,Nfix ($\frac{\mathrm{d}{\mathrm{N}}_{\text{int},\text{Nfix}}}{\mathrm{d}t}$; [$\mathrm{kg}\mathrm{N}{\mathrm{indiv}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{\mathrm{N}}_{\text{int},\text{Nfix}}}{\mathrm{d}t}}=& {\displaystyle }{\mathrm{N}}_{\text{Nfix}}+{\displaystyle \frac{{\mathrm{C}}_{\text{alloc},\text{Nfix}}}{{\mathrm{C}:\mathrm{N}}_{\text{alloc}}}}-{\displaystyle \frac{\mathrm{1}}{{\mathrm{C}:\mathrm{N}}_{\text{Nfix}}}}\\ {\displaystyle }& {\displaystyle }\cdot \left(G_{\text{Nfix}}-{\mathit{\xi }}_{\text{Nfix}}{B}_{\text{Nfix}}-{\text{cost}}_{\text{Nfix}}{\mathrm{N}}_{\text{Nfix}}\right)\\ \text{(A54)}& {\displaystyle }& {\displaystyle }-r_{\text{up},\text{veg}}{\mathrm{N}}_{\text{int},\text{Nfix}}.\end{array}$$

U is calculated as

$$\begin{array}{}\text{(A55)}& U=\sum _k{\mathrm{N}}_{\text{passive}}\left(k\right)+\sum _k{\mathrm{N}}_{\text{active}}\left(k\right)+\sum _j{r}_{\text{up},\text{veg}}{\mathrm{N}}_{\text{int},j}.\end{array}$$

Processes in Sect. A6 occur on the fast timescale (30 min). The purpose of C_int,j and N_int,j is to translate between the fast timescale (30 min) and the daily timescale. Plant C allocation to symbionts occurs on the daily timescale (alongside plant growth), but plant N uptake occurs on the fast timescale.

A8 Dynamic plant C allocation to N uptake relative to plant growth

The order of plant C allocation to growth, symbionts, and rhizosphere priming is determined by C limitation relative to N limitation (Cheng et al., 2014; Finzi et al., 2015; Poorter et al., 2012; Treseder, 2004; Zheng et al., 2019). If a plant is more C-limited than N-limited, $\text{NSC}<\text{NSN}\cdot {\mathrm{C}:\mathrm{N}}_{\text{leaf}}$. The plant allocates C to growth, then to N-fixing bacteria (if associated) and EM (if associated), then to rhizosphere priming, and finally to AM. If a plant is more N-limited than C-limited, $\text{NSC}>\text{NSN}\cdot {\mathrm{C}:\mathrm{N}}_{\text{leaf}}$. The plant allocates C to N-fixing bacteria (if associated) and EM (if associated), then to rhizosphere priming, then to growth, and finally to AM.

A9 Soil N₂O and NO emissions

Soil N₂O and NO emissions occur during nitrification (aerobic oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ with oxygen as an electron acceptor, which produces N₂O and NO as by-products) and denitrification (anaerobic oxidation of organic C with ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as an electron acceptor, which produces N₂O as a by-product).

Nitrification rate in soil layer k (nit(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }\text{nit}\left(k\right)=& {\displaystyle }{\displaystyle \frac{V_{\text{nit},\text{max},\text{ref}}}{\exp \left(\frac{-E_{\mathrm{a},\text{nit}}}{R{T}_{\text{ref}}}\right)}}\exp \left({\displaystyle \frac{-E_{\mathrm{a},\text{nit}}}{RT\left(k\right)}}\right)\\ \text{(A56)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{{\left(\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{3}}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{2.5}}}{f_{\mathit{\theta },\text{max}}}}{\mathrm{NH}}_{\mathrm{4}}\left(k\right),\end{array}$$

where $V_{\text{nit},\text{max},\text{ref}}$ is the maximum nitrification rate, and E_a,nit is the activation energy of nitrification. This follows LM3-SNAP.

Denitrification rate in soil layer k (denit(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{}\text{(A57)}& \text{denit}\left(k\right)=f_{\text{denit}}\sum _i{D}_{\mathrm{C},i,\text{denit}}\left(k\right).\end{array}$$

This follows LM3-SNAP.

Soil N₂O emission rate in soil layer k (N₂O(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\left(k\right)=& {\displaystyle }\mathrm{0.004}\text{nit}\left(k\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{1}/[\mathrm{1}+max\left(\mathrm{3.52},\mathrm{22}\exp \left(-\mathrm{0.8}{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}\left(k\right)}{\text{HR}\left(k\right)}}\right)\right)\\ \text{(A58)}& {\displaystyle }& {\displaystyle }\cdot max\left(\mathrm{0.1},\mathrm{0.015}{\displaystyle \frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}}-\mathrm{0.32}\right)\text{denit}\left(k\right)],\end{array}$$

where and HR(k) is heterotrophic respiration in soil layer k (summation of maintenance respiration, overflow respiration, and decomposition respiration) [$\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]. This follows LM3V-N (Huang and Gerber, 2015).

Soil NO emission rate in soil layer k (NO(k); [$\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{NO}\left(k\right)=\mathrm{0.004}\left(\mathrm{15.2}+{\displaystyle \frac{\mathrm{35.5}{\tan }^{-\mathrm{1}}\left(\mathrm{0.68}\mathit{\pi }\left(\mathrm{2.09}{\left(\mathrm{1}-\frac{\mathit{\theta }\left(k\right)}{{\mathit{\theta }}_{\text{sat}}}\right)}^{\mathrm{4}/\mathrm{3}}-\mathrm{1.68}\right)\right)}{\mathit{\pi }}}\right)\\ \text{(A59)}& {\displaystyle }& {\displaystyle }\cdot \text{nit}\left(k\right).\end{array}$$

This follows LM3V-N (Huang and Gerber, 2015).

Processes in Sect. A8 occur on the fast timescale (30 min).

A10 Additional N sources

N deposition enters N_{U,labile plant-derived}, NH₄, and NO₃. Processes in Sect. A9 occur on the fast timescale (30 min).

Appendix B: Data availability and processing

The US Forest Inventory and Analysis (FIA) database (US Forest Service, 2020a) was downloaded in 2019 from https://apps.fs.usda.gov/fia/datamart/CSV/datamart_csv.html, last access: 7 November 2019. For tree data, plots with Plot Design Code 1 (national plot design) were selected, dead and cut trees were excluded, trees measured at the root collar were excluded, trees with visually estimated or modelled heights were excluded, and accessible forest land was selected (excludes agriculture and urban areas). Canopy trees had crown class code open grown, dominant, or codominant. Understory trees had crown class code intermediate or overtopped. For soil data, mineral soil layers were selected. For seedling data, plots with Plot Design Code 1 (national plot design) were selected, and accessible forest land was selected (excludes agriculture and urban areas).

The US FIA Forest Health Monitoring database (US Forest Service, 2020b) was downloaded in 2019 from https://www.fia.fs.fed.us/tools-data/other_data/index.php, last access: November 2019. Dead trees and trees with damaged crowns were excluded. Timberland and woodland land uses were selected.

The Biomass and Allometry Database (BAAD; Falster et al., 2015) was downloaded from https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/14-1889.1, last access: November 2019. Field/wild plants were selected. Temperate forest plants were selected.

TRY database data (Kattge et al., 2020) were downloaded in 2019 from https://www.try-db.org/TryWeb/Prop0.php, last access: November 2019.

Total soil C, total soil N, soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$, soil nitrate ${\mathrm{NO}}_{\mathrm{3}}^{-}$, N mineralization rate, net nitrification rate, and bulk density from CHL (Knoepp, 2009a, b, 2018) were downloaded from https://coweeta.uga.edu/ResearchData.html, last access: November 2019.

Soil N₂O and NO emissions are from Stehfest and Bouwman (2006). Data from deciduous temperate forests were selected.

Net ecosystem production at the hourly timescale (Oishi, 2020) were downloaded from https://ameriflux.lbl.gov/doi/AmeriFlux/US-Cwt, last access: November 2019. Growing season was May through October.

Gross primary production, net primary production, heterotrophic respiration and net ecosystem production are from the ForC database (Anderson-Teixeira et al., 2018). Data from deciduous temperate broadleaf forests were selected.

Appendix C: Model parameterization

Equation (A40) was fit using nonlinear least squares with the US FIA Forest Health Monitoring database (US Forest Service, 2020b) for each species to calculate α_CA and θ_CA (Figs. D10 and D11). Equation (A41) was fit using nonlinear least squares with the US FIA database tree data (US Forest Service, 2020a) for each species to calculate α_HT, θ_HT, and γ_HT (Figs. D12 and D13). Equation (A42) was fit using nonlinear least squares with the US FIA database tree data (US Forest Service, 2020a) for each species to calculate α_BM (Figs. D14 and 15). Note that wood mass (aboveground and belowground) is calculated by the US FIA using allometric equations of dbh.

The following equations from Weng et al. (2015) were utilized with the US FIA database tree data from consecutive censuses for each species (US Forest Service, 2020a) to calculate canopy tree background mortality rate (μ_C) and understory tree background mortality rate (μ_U):

$$\begin{array}{}\text{(C1)}& {\displaystyle }& {\displaystyle }{n}_{\text{dead},\mathrm{C}}=n_{\mathrm{C}}\left(\mathrm{1}-e^{-{\mathit{\mu }}_{\mathrm{C}}}\right),\text{(C2)}& {\displaystyle }& {\displaystyle }{n}_{\text{dead},\mathrm{U}}=n_{\mathrm{U}}\left(\mathrm{1}-e^{-{\mathit{\mu }}_{\mathrm{U}}\left(\mathrm{1}+\frac{\mathrm{4}{e}^{\mathrm{30}\left(\mathrm{0.025}-D_{\mathrm{U}}\right)}}{\mathrm{1}+e^{\mathrm{30}\left(\mathrm{0.025}-D_{\mathrm{U}}\right)}}\right)}\right),\end{array}$$

where n_dead,C is the mean number of dead canopy individuals in a given year [indiv m⁻²], n_C is the mean number of canopy individuals in a given year [indiv m⁻²], n_dead,U is the mean number of dead understory individuals in a given year [indiv m⁻²], n_U is the mean number of understory individuals in a given year [indiv m⁻²], and D_U is the mean diameter at breast height of an understory individual in a given year [m].

The following equation from Martinez Cano et al. (2020) was fit using nonlinear least squares with BAAD (Falster et al., 2015) to calculate the proportionality constant between crown area and cross-sectional sapwood area (φ_CSA) (Figs. D16 and D17):

$$\begin{array}{}\text{(C3)}& {\mathit{\phi }}_{\text{CSA}}={\displaystyle \frac{S{A}_{\text{sw}}\left({\mathit{\gamma }}_{\text{HT}}+D^{{\mathit{\theta }}_{\text{HT}}}\right)}{D^{{\mathit{\theta }}_{\text{CA}}+{\mathit{\theta }}_{\text{HT}}}}},\end{array}$$

where SA_sw is the cross-sectional sapwood area at breast height [m²], and D is the diameter [m] from BAAD (for deciduous angiosperm).

The following equation from Weng et al. (2015) was utilized with BAAD (Falster et al., 2015) to calculate the proportionality constant between root area and leaf area (φ_RL):

$$\begin{array}{ll}{\displaystyle }{\mathit{\phi }}_{\text{RL}}& {\displaystyle }={\displaystyle \frac{{\mathrm{C}}_{\mathrm{r}}}{{\mathrm{C}}_{\mathrm{l}}}}\cdot \text{LMA}\cdot \text{SRA}\\ {\displaystyle }& {\displaystyle }=\mathrm{0.43}\cdot \mathrm{0.0375}\left[\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}\right]\mathrm{45}\left[{\mathrm{m}}^{\mathrm{2}}(\mathrm{kg}\mathrm{C}{)}^{-\mathrm{1}}\right]=\mathrm{0.79},\end{array}$$

where 0.43 is from BAAD (for deciduous angiosperm), 0.0375 kg C m⁻² is from Poorter et al. (2009) (for deciduous trees), and 45 m² kg C is from Jackson et al. (1997) (for temperate deciduous forests).

The following equation from Martinez Cano et al. (2020) was utilized with BAAD (Falster et al., 2015) to calculate the fraction of sapwood in branches (f_br)

$$f_{\text{br}}={\displaystyle \frac{{\mathrm{C}}_{\text{br}}}{{\mathrm{C}}_{\text{sw}}}}=\mathrm{0.255},$$

where 0.255 is from BAAD. Note that we did not distinguish between deciduous angiosperm and evergreen gymnosperm due to insufficient data.

The following equation (from Eq. A20) was used to calculate r_{Nfix asymb}:

$$\begin{array}{ll}{\displaystyle }{r}_{\text{Nfix asymb}}& {\displaystyle }={\displaystyle \frac{{\mathrm{N}}_{\text{Nfix asymb}}}{{\mathrm{C}}_{\mathrm{M}}}}={\displaystyle \frac{\mathrm{12}\left[\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}\right]}{\mathrm{50}\mathrm{000}\left[\mathrm{kg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}\right]\times \mathrm{0.01}}}\\ {\displaystyle }& {\displaystyle }=\mathrm{0.024}\left[\mathrm{kg}\mathrm{N}(\mathrm{kg}\mathrm{C}\mathrm{yr}{)}^{-\mathrm{1}}\right],\end{array}$$

where 12 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ is from Reed et al. (2011), 50 000 kg C ha⁻¹ is from Scharlemann et al. (2014), and 1 % (0.01) is from Chapin et al. (2011).

The following equation (from Eq. A38) was used to calculate r_leakage,C:

$$\begin{array}{ll}{\displaystyle }{r}_{\text{leakage},\mathrm{C}}& {\displaystyle }={\displaystyle \frac{L_{\mathrm{C}}}{\text{NSC}}}={\displaystyle \frac{\mathrm{0.05}\cdot \text{PP}}{q\cdot {\mathrm{C}}_{\mathrm{l}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{0.05}\cdot \mathrm{13}\mathrm{050}\left[\mathrm{kg}\mathrm{C}(\mathrm{ha}\mathrm{yr}{)}^{-\mathrm{1}}\right]}{\mathrm{4}\cdot \mathrm{3.2}\left[\mathrm{kg}\mathrm{C}{\mathrm{indiv}}^{-\mathrm{1}}\right]\cdot \mathrm{390}\left[\mathrm{indiv}{\mathrm{ha}}^{-\mathrm{1}}\right]}}\\ {\displaystyle }& {\displaystyle }=\mathrm{0.131}\left[{\mathrm{yr}}^{-\mathrm{1}}\right],\end{array}$$

where 0.05 is from Jones et al. (2009), 13 050 $\mathrm{kg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ is from the ForC database (Anderson-Teixeira et al., 2018), 3.2 kg C indiv⁻¹ is from BAAD (Falster et al., 2015), and ∼390 indiv ha⁻¹ is from Crowther et al. (2015) (for temperate forests). Note that we did not distinguish between deciduous angiosperm and evergreen gymnosperm in BAAD because other data are from all temperate forests.

r_growth from LM3-SNAP was adapted to the units of LM4.1:

$$\begin{array}{ll}{\displaystyle }{r}_{\text{growth}}& {\displaystyle }=\mathrm{0.3}\left[\mathrm{kg}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}\right]\left[\mathrm{10}\mathrm{000}{\mathrm{m}}^{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}\right]\\ {\displaystyle }& {\displaystyle }/(\mathrm{8}\left[\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}\right]\left[\mathrm{kg}\mathrm{C}(\mathrm{2}\mathrm{kg}{)}^{-\mathrm{1}}\right],\\ {\displaystyle }& {\displaystyle }\mathrm{106}\left[\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}\right]\left[\mathrm{kg}\mathrm{C}(\mathrm{2}\mathrm{kg}{)}^{-\mathrm{1}}\right],\\ {\displaystyle }& {\displaystyle }\mathrm{4}\left[\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}\right]\left[\mathrm{kg}\mathrm{C}(\mathrm{2}\mathrm{kg}{)}^{-\mathrm{1}}\right],\\ {\displaystyle }& {\displaystyle }\mathrm{250}\left[\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}\right]\left[\mathrm{kg}\mathrm{C}(\mathrm{2}\mathrm{kg}{)}^{-\mathrm{1}}\right])\\ {\displaystyle }& {\displaystyle }=\mathrm{65}{\mathrm{yr}}^{-\mathrm{1}},\end{array}$$

where 0.3 $\mathrm{kg}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$ is from LM3-SNAP; 8, 106, and 4 kg ha⁻¹ are from Boring and Swank (1984); and 250 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ is from Binkley et al. (1992).

r_Nfix was calculated with data from Bytnerowicz et al. (2021), which measured nodule biomass % C and symbiotic BNF rate per unit nodule biomass using ¹⁵N₂ incorporation.

Parameters related to plant C allocation to symbionts (f_alloc,AM, f_alloc,EM, and f_alloc,Nfix) were estimated by examining a range of values (0.01, 0.02, 0.05, 0.10, 0.20, 0.50, and 1.00). Symbiont biomass C was compared to literature observations (that were not used for model parameterization) from Boring and Swank (1984) for nodule biomass C and Zhu and Miller (2003) for mycorrhizal biomass C. Parameter values were chosen such that simulated symbiont biomass C was most similar to observed symbiont biomass C.

Appendix D: Appendix D figures and tables

Figure D1 Simulated dbh growth rate for (a) Acer/non-fixer and (b) Robinia compared to FIA data (in North Carolina). Simulated data are trees with dbh>12.7 cm to reflect the dbh range of FIA data and are weighted by the stand age distribution of FIA data (Fig. D18). Simulated and FIA data are scaled to display an equal maximum density. LM3-SNAP cannot distinguish between plant cohorts with different dbh's or distinguish between Robinia and Acer.

Figure D2 Simulated dbh distribution compared to FIA data (in North Carolina). Simulated data are trees with dbh>12.7 cm to reflect the dbh range of FIA data and are weighted by the stand age distribution of FIA data (Fig. D18). LM3-SNAP cannot distinguish between plant cohorts with different dbh's or distinguish between Robinia and Acer.

Figure D3 Controls of asymbiotic BNF. (a) Soil microbe biomass C, (b) soil microbe biomass N, (c) total soil C, and (d) total soil N.

Figure D4 Controls of symbiotic BNF. (a) Nodule biomass C, (b) non-structural C, and (c) N stress.

Figure D5 Simulated soil C and N pools, soil N fluxes, and soil N loss rates from LM4.1-BNF, LM4.1-BNFNPP, and LM4.1-BNFET. (a) Simulated total soil C (depth 0–10 cm). (b) Simulated total soil N (depth 0–10 cm). (c) Simulated soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (depth 0–10 cm). (d) N mineralization rate and net nitrification rate (depth 0–10 cm). (e) Simulated N₂O and NO emission rate and simulated dissolved organic N (DON), ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the site data which are from mature forests. Error bars indicate 2 standard deviations.

Figure D6 Simulated gross primary production (GPP), net primary production (NPP), heterotrophic respiration (HR), and net ecosystem production (NEP) from LM4.1-BNF, LM4.1-BNFNPP, and LM4.1-BNFET. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the data which are from mature forests. Error bars indicate 2 standard deviations.

Figure D7 Simulated soil C and N pools, soil N fluxes, and soil N loss rates from LM4.1-BNF initialized with both Robinia and Acer, only Acer, and only N-fixer Acer, with and without asymbiotic BNF. (a) Simulated total soil C (depth 0–10 cm). (b) Simulated total soil N (depth 0–10 cm). (c) Simulated soil ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (depth 0–10 cm). (d) N mineralization rate and net nitrification rate (depth 0–10 cm). (e) Simulated N₂O and NO emission rate and simulated dissolved organic N (DON), ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rate. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the site data which are from mature forests. Error bars indicate 2 standard deviations.

Figure D8 Simulated gross primary production (GPP), net primary production (NPP), heterotrophic respiration (HR), and net ecosystem production (NEP) from LM4.1-BNF initialized with both Robinia and Acer, only Acer, and only N-fixer Acer, with and without asymbiotic BNF. Simulated data are averaged over the last 100 of the 300 years of simulation to reflect the data which are from mature forests. Error bars indicate 2 standard deviations.

Figure D9 Simulated absolute basal area of Acer and Robinia over time compared to FIA data (in North Carolina). Simulated data are trees with dbh>12.7 cm to reflect the dbh range of FIA data. FIA data of all non-fixing trees are aggregated to represent Acer. Each point represents an FIA plot.

Figure D10 Crown area model fit to FIA Forest Health Monitoring data for Acer.

Figure D11 Crown area model fit to FIA Forest Health Monitoring data for Robinia.

Figure D12 Height model fit to FIA data for Acer.

Figure D13 Height model fit to FIA data for Robinia.

Figure D14 Wood mass model fit to FIA data for Acer (displayed at mean height of Acer).

Figure D15 Wood mass model fit to FIA data for Robinia (displayed at mean height of Robinia).

Figure D16 Sapwood area at breast height model fit to BAAD for Acer.

Figure D17 Sapwood area at breast height model fit to BAAD for Robinia.

Figure D18 Histogram of FIA stand age for FIA plots in North Carolina. These distributions were used to weigh the data displayed in Figs. 2 and 3.

Table D1 Vegetation-type-specific parameters.

Table D2 General parameters.

Table D3 Summary of the spin-up and numerical experiments.

Table D4 Initial densities and heights for simulations. The initial densities of Robinia and Acer were derived from the US FIA database seedling data for plots with at least one Robinia individual in North Carolina. Diameter at breast height is determined from height by allometry (Eq. A41).

Table D5 Validated variables and data sources.