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 kgNha−1 yr−1) in comparison to observations (11 kgNha−1 yr−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 kgNha−1 yr−1 at 15 years and then declines to ∼ 0 kgNha−1 yr−1 at 300 years, similarly to observed symbiotic BNF, which reaches 75 kgNha−1 yr−1 at 17 years and then declines to ∼ 0 kgNha−1 yr−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 CO2 concentration and other global change factors.


4144
S. Kou-Giesbrecht et al.: A novel representation of biological nitrogen fixation plant growth will constrain the terrestrial C sink under elevated atmospheric CO 2 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;, 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 downregulated under non-N limitation of plant growth . BNF could, as such, be pivotal to overcoming N limitation of plant growth under elevated atmospheric CO 2 concentration (Liang et al., 2016;Terrer et al., 2016Terrer et al., , 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 metaanalysis 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 , represent a single general plant C pool capable of BNF and cannot represent community dynamics. In observed ecosystems, sym-biotic 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 downregulate 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 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 Nfixing 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.
4. N limitation. N limitation is determined by current stored non-structural N relative to the demand for nonstructural 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

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) andMartinez Cano et al. (2020) and are not directly influenced by N availability but are indirectly influenced by N availability via its effect on growth.  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 (NH + 4 ) and nitrate (NO − 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.
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 vegetationtype-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.

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  (Eq. A8).
The symbiotic BNF rate by N-fixing bacteria (N Nfix ; 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). Respiration associated with symbiotic BNF is cost Nfix N 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.

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 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 −2 ], and f a (T (k)) is the soil temperature dependence function.
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 −1 ]) 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 ( dNSC dt ; [kg C indiv −1 yr −1 ]) is where P is the photosynthesis rate [kg C indiv −1 yr −1 ], R is the respiration rate (maintenance and growth) [kg C indiv −1 yr −1 ], G C,l is the growth rate of the leaf C pool (C l ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], G C,r is the growth rate of the fine-root C pool (C r ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], C sw is the growth rate of the sapwood C pool (C sw ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], G C,seed is growth rate of the seed C pool (C seed ; [kg C indiv −1 ]) [kg C indiv −1 yr −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 −1 ]) 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 ( dNSN dt ; [kg N indiv −1 yr −1 ]) is where U is the N uptake rate via roots and symbionts (Eq. A55) [kg N indiv −1 yr −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) andMartinez Cano et al. (2020). The total allocation of NSC to growth is determined by the target NSC (NSC target ; [kg C indiv −1 ]) 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 −1 ]), which reflects the ability of a plant to refoliate after defoliation (Hoch et al., 2003;Richardson et al., 2013), and is calculated as 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 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 [kg C indiv −1 yr −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 Nlimited (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 where NSN target is the target NSN [kg N indiv −1 ]. 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 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 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 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); [kg N indiv −1 yr −1 ]) increases with N stress and is calculated as where f (NO 3 (k), NH 4 (k), f rhiz (k)) is a function of the NH + 4 pool in soil layer k (NH 4 (k); [kg N m −2 ]), the NO − 3 pool in soil layer k (NO 3 (k); [kg N m −2 ]), and the rhizosphere volume fraction of soil layer k (f rhiz (k); [m 3 m −3 ]) and is given in Eq. (A3).

Plant C allocation to symbionts (AM, EM, and N-fixing bacteria)
The rate of C allocation to AM (C alloc,AM ; 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 . The rate of C allocation to EM (C alloc,EM ; [kg C indiv −1 yr −1 ]) is 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 . 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 ; [kg C indiv −1 yr −1 ]) is where f alloc,Nfix is the fraction of NSC allocated to Nfixing bacteria per unit time, and f alloc,Nfix,min is the minimum fraction of NSC allocated to N-fixing bacteria per unit time. For plants with a perfectly facultative BNF strategy f alloc,Nfix,min = 0, and for plants with an incomplete downregulator BNF strategy f alloc,Nfix,min > 0. Robinia has an incomplete down-regulator BNF strategy. For plants with an obligate BNF strategy, C alloc,Nfix is 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; [kg N indiv −1 yr −1 ]) is 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).

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, NSC < NSN·C : N 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 Climited, NSC > NSN · C : N leaf , the plant allocates C to Nfixing bacteria (if associated) and EM (if associated), then to rhizosphere priming, then to growth, and finally to AM.

Soil N 2 O and NO emissions
Soil N 2 O and NO emissions occur during nitrification (aerobic oxidation of NH + 4 with oxygen as an electron acceptor, which produces N 2 O and NO as by-products) and denitrification (anaerobic oxidation of organic C with NO − 3 as an electron acceptor, which produces N 2 O as a by-product).
Following LM3V-N (Huang and Gerber, 2015), , θ sat is saturation volumetric soil water content, and denit(k) is denitrification rate in soil layer k [kg N m −2 yr −1 ] (Eq. A57). Following LM3V-N (Huang and Gerber, 2015), soil NO emission rate in soil layer k where NPP is net primary production [kg C m −2 yr −1 ], and a NPP and b NPP are constants from Meyerholt et al. (2016).
where ET is evapotranspiration [mm yr −1 ], and a ET and b ET are constants from Meyerholt et al. (2016). BNF NPP and BNF ET enter NH 4 (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.

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 preindustrial atmospheric CO 2 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 spinup) and the LM4.1-BNF spin-up soil C and N pools. We ran numerical experiments for another 300 years at current atmospheric CO 2 concentration (324.53 ppm). See Table D3 for a summary of atmospheric CO 2 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.
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 preindustrial atmospheric CO 2 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 2 concentration. See Table D3 for a summary of atmospheric CO 2 concentration, meteorological forcings, and N deposition rates used in the spin-up and numerical experiments. We initialized the LM3-SNAP spin-up and the 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).
LM3-SNAP numerical experiments with a temperate deciduous vegetation type.

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).
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 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. ( Fig. 6; US Forest Service, 2020a), site observations from CHL ( Fig. 6; Knoepp, 2009aKnoepp, , b, 2018Swank 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  Table D5 for a summary of the validated variables and data sources.

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 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 NH + 4 and NO − 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 2 O and NO emission rates compared to a meta-analysis estimate for temperate forests and simulated dissolved organic N (DON), NH + 4 , and NO − 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 2 O or NO emissions. 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 kg N ha −1 yr −1 over the final 100 years), which is not represented in LM3-SNAP (Fig. 3a), in comparison to observations from CHL (11 kg N ha −1 yr −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 kg N ha −1 yr −1 at 15 years and then declines to ∼ 0 kg N ha −1 yr −1 at 300 years (Fig. 3b). LM3-SNAP simulates high symbiotic BNF rate in late succession (mean 8 kg N ha −1 yr −1 over the final 100 years; Fig. 3b). Observations from CHL suggest that symbiotic BNF which reaches 75 kg N ha −1 yr −1 at 17 years then declines to ∼ 0 kg N ha −1 yr −1 in late suc-cessional 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).
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 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. 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 Mg C ha −1 yr −1 respectively). LM4.1-BNF overestimates HR (mean 7.5 vs. 4.7 Mg C ha −1 yr −1 ) and consequentially underestimates net ecosystem production (NEP) (mean 0.6 vs. 4.8 Mg C ha −1 yr −1 ) in comparison to observations.

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, Nfixing 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 −1 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 NO − 3 leaching rate (mean 8.3 vs. 0.1 kg N ha −1 yr −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).

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 −1 over the final 100 years), followed by LM4.1-BNF ET (mean 70 Mg C ha −1 over the final 100 years) and LM4.1-BNF NPP (mean 0 Mg C ha −1 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 kg N ha −1 yr −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 kg N ha −1 yr −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 pi- oneer 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. 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.

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 ini- tialized with both Robinia and Acer seedlings (mean 48.6 vs. 80.9 Mg C ha −1 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 kg N ha −1 yr −1 ), reaching a similar level to LM4.1-BNF simulations initialized with both Robinia and Acer seedlings after 100 years. 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 −1 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 kg N ha −1 yr −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 Nfixer 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).

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 −1 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 kg N ha −1 yr −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.

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).

S. Kou-Giesbrecht et al.: A novel representation of biological nitrogen fixation
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.

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 restructuring 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.

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 Nfixing 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 2 O emissions in LM4.1-BNF enables the estimation of the magnitude of the terrestrial NO and N 2 O source, which can be driven by BNF (Kou-Giesbrecht and . 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 2 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 2 sequestration by terrestrial ecosystems under elevated atmospheric CO 2 concentration.

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

A2.1 Active root uptake of inorganic N
The rate of active root uptake of inorganic N in soil layer k where r NO 3 and r NH 4 are rate constants, f rhiz (k) is the rhizosphere volume fraction of soil layer k [m 3 m −3 ], z(k) is the thickness of soil layer k, k M,NO 3 and k M,NH 4 are half-saturation constants, and n indiv is the spatial density [indiv m −2 ]. f rhiz (k) is calculated as 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 −1 ], 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 NO − 3 and NH + 4 uptake by AM in soil layer k (N AM,NO 3 (k) and N AM,NH 4 (k) respectively; where r NO 3 ,AM and r NH 4 ,AM are rate constants, k AM,NO 3 and k AM,NH 4 are half-saturation constants, B AM (k) is the biomass C of AM in soil layer k [kg C indiv −1 ], and k AM is a halfsaturation 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 C,i,EM (k) and D N,i,EM (k) respectively; [kg C indiv −1 yr −1 ] and [kg N indiv −1 yr −1 ]), are The rates of C and N uptake by EM in soil layer k (C EM (k) and N EM (k) respectively; [kg C indiv −1 yr −1 ] and [kg N indiv −1 yr −1 ]) are where ε C,i,EM is the C uptake efficiency of soil C type i by EM, ε N,i,EM is the N uptake efficiency of soil N type i by EM, r NO 3 ,EM and r NH 4 ,EM are rate constants, k EM,NO 3 and k EM,NH 4 are half-saturation constants, and k EM is a half- . 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 ; [kg N indiv −1 yr −1 ]) is 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 plantderived, labile microbe-derived, or recalcitrant, (D C,i (k) and D N,i (k) respectively; [kg C m −2 yr −1 ] and [kg N m −2 yr −1 ]), are where V max,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 −2 ], 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 C,i,denit,pot (k) and D N,i,denit,pot (k) respectively; [kg C m −2 yr −1 ] and [kg N m −2 yr −1 ]) are where V denit,max,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 C,i,denit (k) and D N,i,denit (k) respectively; [kg C m −2 yr −1 ] and [kg N m −2 yr −1 ]) are where k denit is the half-saturation constant, and f denit is the stoichiometric ratio of NO − 3 demand for C decomposition. D C,i,denit (k) and D N,i,denit (k) include NO − 3 limitation of denitrification. This follows LM3-SNAP.
The rate of change of biomass C of soil microbes in soil 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 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); [kg C m −2 yr −1 ]) is This is released as CO 2 . This follows LM3-SNAP. The asymbiotic BNF rate of soil microbes in soil layer k where r Nfix asymb is a rate constant and f a (T (k)) is the soil temperature dependence function. 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
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.
. 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 (φ N (k) ≥ − Imm max (k)) or N-limited (φ N (k) < − Imm max (k)) and are calculated as follows.
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); [kg N m −2 yr −1 ]) is where V max,ref,NH 4 is the maximum NH + 4 immobilization rate, E a,NH 4 is the activation energy of NH + 4 immobilization, V max,ref,NO 3 is the maximum NO − 3 immobilization rate, and S. Kou-Giesbrecht et al.: A novel representation of biological nitrogen fixation E a,NO 3 is the activation energy of NO − 3 immobilization. This follows LM3-SNAP.
When growth of soil microbes in soil layer k is N-limited (φ N (k) < − Imm max (k)), there is overflow respiration of excess C. The overflow respiration of excess C in soil layer k (R overflow (k); [kg C m −2 yr −1 ]) is This is released as CO 2 . If φ N (k) > 0, the net N mineralization flux in soil layer k is The fraction of the net N immobilization flux in soil layer k that is NH + 4 immobilization is .

The fraction of the net N immobilization flux in soil layer
. 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 −1 ]) 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 ((dNSC)/(dt); [kg C indiv −1 yr −1 ]) is where P is the photosynthesis rate [kg C indiv −1 yr −1 ], R is the respiration rate (maintenance and growth) [kg C indiv −1 yr −1 ], G C,l is the growth rate of the leaf C pool (C l ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], G C,r is the growth rate of the fine-root C pool (C r ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], G C,sw is the growth rate of the sapwood C pool (C sw ; [kg C indiv −1 ]) [kg C indiv −1 yr −1 ], G C,seed is the growth rate of the seed C pool (C seed ; [kg C indiv −1 ]) [kg C indiv −1 yr −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 −1 ]) 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 ((dNSN)/(dt); [kg N indiv −1 yr −1 ]) is where U is the N uptake rate via roots and symbionts (Eq. A55) [kg N indiv −1 yr −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 −1 ]) 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 −1 ]), which reflects the ability of a plant to refoliate after defoliation (Hoch et al., 2003;Richardson et al., 2013), and is calculated as 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 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 [kg C indiv −1 yr −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 where NSN target is the target NSN [kg N indiv −1 ]. 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 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 −1 ]) is If N demand > 0.5NSN, G C,l , G C,r , G C,sw , and G C,seed are reduced to prevent N demand from depleting NSN 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 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 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.
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) where α HT , θ HT , and γ HT are allometry parameters. Wood mass (W ; [kg C]) is a function of diameter at breast height and height: 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 ; 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 . The rate of C allocation to EM (C alloc,EM ; [kg C indiv −1 yr −1 ]) is C alloc,EM = f alloc,EM NSC N stress , 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 . 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 ; [kg C indiv −1 yr −1 ]) is where f alloc,Nfix is the fraction of NSC allocated to Nfixing bacteria per unit time and f alloc,Nfix,min is the minimum fraction of NSC allocated to N-fixing bacteria per unit time. For plants with a perfectly facultative BNF strategy f alloc,Nfix,min = 0, and for plants with an incomplete downregulator BNF strategy f alloc,Nfix,min > 0. We model Robinia with an incomplete down-regulator BNF strategy.
For plants with an obligate BNF strategy, C alloc,Nfix is C alloc,Nfix = f alloc,Nfix NSC.
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; [kg N indiv −1 yr −1 ]) is 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 −1 ]). The rate of change of C int,j ((dC int,j )/(dt); 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; [kg C indiv −1 yr −1 ]) is where r growth is the growth rate of a symbiont. The growth rate of EM (G EM ; [kg C indiv −1 yr −1 ]) is EM, or N-fixing bacteria) is released as CO 2 (growth respiration). The rate of change of biomass C of a symbiont with strategy j ((dB j )/(dt); j = AM or EM; [kg C indiv −1 yr −1 ]) is 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 . ξ j B j is released as CO 2 (maintenance respiration). (B j )/(τ j ) enters C U,labile microbe-derived and (B j )/(C : N j τ j ) 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 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. ξ Nfix B Nfix is released as CO 2 (maintenance respiration), and cost Nfix N Nfix is released as CO 2 (respiration associated with symbiotic BNF). (B Nfix )/(τ Nfix ) enters C U,labile microbe-derived , and (B Nfix )/(C : N Nfix τ Nfix ) enters N U,labile microbe-derived , where C : N Nfix is the C : N ratio of Nfixing bacteria. N acquired by symbionts is transferred to an intermediate N pool (N int,j ; [kg N indiv −1 ]). The rate of change of N int,j ((dN int,j )/(dt); j = AM or EM; [kg N indiv −1 yr −1 ]) is where r up,veg is the rate of plant N uptake from N int,j . The rate of change of N int,Nfix ( U is calculated as 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, NSC < NSN · C : N 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 Climited, NSC > NSN · C : N leaf . The plant allocates C to Nfixing bacteria (if associated) and EM (if associated), then to rhizosphere priming, then to growth, and finally to AM.

A9 Soil N 2 O and NO emissions
Soil N 2 O and NO emissions occur during nitrification (aerobic oxidation of NH + 4 with oxygen as an electron acceptor, which produces N 2 O and NO as by-products) and denitrification (anaerobic oxidation of organic C with NO − 3 as an electron acceptor, which produces N 2 O as a by-product).
Nitrification rate in soil layer k (nit(k); [kg N m −2 yr −1 ]) is where V nit,max,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); where and HR(k) is heterotrophic respiration in soil layer k (summation of maintenance respiration, overflow respiration, and decomposition respiration) [kg C m −2 yr −1 ]. This follows LM3V-N (Huang and Gerber, 2015). Soil NO emission rate in soil layer k (NO(k); This follows LM3V-N (Huang and Gerber, 2015 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 ):  (Falster et al., 2015) to calculate the proportionality constant between crown area and cross-sectional sapwood area (ϕ CSA ) (Figs. D16 and D17): where SA sw is the cross-sectional sapwood area at breast height [m 2 ], 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 ):  (Falster et al., 2015) to calculate the fraction of sapwood in branches (f br ) 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 : where 12 kg N ha −1 yr −1 is from Reed et al. (2011), 50 000 kg C ha −1 is from Scharlemann et al. (2014), and 1 % (0.01) is from Chapin et al. (2011).
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 15 N 2 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.  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 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.

BNF representation Initialized species
Data availability. Data sets are described and listed in Appendix B.
Author contributions. SKG designed and implemented LM4.1-BNF in collaboration with DNLM, SM, and ES. IMC, SWP, and TAB advised on the new model. SM and ES developed the initial implementation of LM4.1 with LM3-SNAP-based N cycling. SKG conducted the analysis and wrote the initial draft. All authors contributed to the writing.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.