Modelled forest ecosystem carbon–nitrogen dynamics with integrated mycorrhizal processes under elevated CO₂

Almost 95 % of all terrestrial plant species form symbioses with mycorrhizal fungi that mediate plant–soil interactions: mycorrhizae facilitate plant nitrogen (N) acquisition and are, therefore, vital for plant growth, but they also build a pathway for plant-assimilated carbon (C) into the rhizosphere. Therefore, mycorrhizae likely play an important role in shaping the response of ecosystems to environmental changes such as rising atmospheric carbon dioxide (CO₂) concentrations, which can increase plant N demand and the transfer of plant C assimilation to the soil. While the importance of mycorrhizal fungi is widely recognised, they are rarely represented in current terrestrial biosphere models (TBMs) explicitly. Here, we present a novel, dynamic plant–mycorrhiza–soil model as part of the QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) TBM. This new model is based on mycorrhizal functional types that either actively mine soil organic matter (SOM) for N or enhance soil microbial activity through increased transfer of labile C into the rhizosphere, thereby (passively) priming SOM decomposition. Using the Duke Free-Air CO₂ Enrichment (FACE) experiment, we show that mycorrhizal fungi can have important effects on projected SOM turnover and plant nutrition under ambient as well as elevated-CO₂ treatments. Specifically, we find that including enhanced active mining of SOM for N in the model allows one to more closely match the observations with respect to observed decadal responses of plant growth, plant N acquisition and soil C dynamics to elevated CO₂, whereas a simple enhancement of SOM turnover by increased below-ground C transfer of mycorrhizae is unable to replicate the observed responses. We provide an extensive parameter uncertainty study to investigate the robustness of our findings with respect to model parameters that cannot readily be constrained by observations. Our study points to the importance of implementing mycorrhizal functionalities in TBMs as well as to further observational needs to better constrain mycorrhizal models and to close the existing major knowledge gaps in actual mycorrhizal functioning.

1 Introduction

The land biosphere sequesters about a quarter of the annual anthropogenic CO₂ emissions, thereby slowing human-induced climate change (Le Quéré et al., 2018). However, it is unclear if this buffering effect will continue into the future (Ciais et al., 2014). Multiple lines of evidence (Norby et al., 2005; Jiang et al., 2020; Walker et al., 2021) and the majority of nitrogen-enabled terrestrial biosphere models (TBMs; Zaehle et al., 2014; Wieder et al., 2015; Fleischer et al., 2019) suggest that CO₂-fertilised plant growth, which is one of the major causes of this uptake, will be reduced by increasing nutrient limitation. Whole-ecosystem experiments, in which plants are exposed to elevated levels of atmospheric CO₂ (eCO₂), show a range of responses from ongoing CO₂-fertilised growth (McCarthy et al., 2007) to progressive nutrient limitation (Norby et al., 2010). One possible explanation for these divergent responses is the differential ability of plants to accesses soil nutrients through mycorrhizal fungi (Terrer et al., 2016).

Mycorrhizal fungi are ubiquitous and form symbioses with almost 95 % of all terrestrial plant species (Read, 1991). It is commonly assumed that this mutualistic relationship between plants and fungi is an evolutionary adaptation to nutrient-poor conditions (Read, 1991), whereby host plants benefit from mycorrhizal nutrient supply (Ames et al., 1983; Marschner and Dell, 1994; Ek et al., 1997). In return, plants export carbohydrates to their symbiotic fungi, covering most of the fungal carbon (C) demand (Gorka et al., 2019). Plant–mycorrhiza symbioses are important for soil C sequestration and, thus, soil organic matter (SOM) content, either via the allocation of C into below-ground tissues such as roots or mycorrhizal biomass or by rhizodeposition, i.e. C exudation into the rhizosphere of the soil (Godbold et al., 2006; Frey, 2019; Gorka et al., 2019).

Mycorrhizae are generally classified based on their morphology into arbuscular mycorrhizae, ectomycorrhizae, ericoid mycorrhizae and orchid mycorrhizae, but the first two are the most common forms (Read, 1991; Johnson and Jansa, 2017). Arbuscular mycorrhizae (AMs) directly enter plant root cells and build symbiotic associations with more than 80 % of all terrestrial plant families, especially with crop and grass species, but also with most tropical trees, which is why they are often assumed to be more important for plant phosphorus (P) nutrition (Read, 1991; Hodge et al., 2001). Ectomycorrhizae (ECMs) are hosted by only 10 % of plant species, including most temperate and boreal tree species and, contrary to previous findings, some tropical tree species (Corrales et al., 2018). This results in ECM interactions covering a large fraction (26.5 %) of mycorrhizal relationships in terrestrial plant biomass (Soudzilovskaia et al., 2019). ECMs are usually regarded as nitrogen (N) nutrition supporters and only cover plant roots, working as a storage and exchange area for carbohydrates and nutrients as well as protection against infections and a physical barrier between the roots and soil nutrients (Read, 1991; Jordy et al., 1998; Laczko et al., 2004; Pritsch et al., 2004; Finlay, 2008).

Given the vast diversity of mycorrhiza species, even after decades of research, there is still no complete understanding on how plant–mycorrhiza symbioses affect ecosystem nutrition or the stabilisation and turnover of SOM (Jansa and Treseder, 2017; Tedersoo and Bahram, 2019). Some species accelerate SOM decay by C exudation to soil, which activates the microbial community and, thus, enhances SOM decomposition rates (Hodge et al., 2001; Cheng et al., 2012; Lindahl and Tunlid, 2015; Paterson et al., 2016; Lang et al., 2021). Other mycorrhizal species have the ability to take up nutrients directly from SOM in addition to mineral nutrients (Wu et al., 2005; Pérez-Tienda et al., 2012; Hodge and Storer, 2015) or even act saprotrophically, i.e. they are able to gain C from SOM (Treseder et al., 2007; Malcolm et al., 2008; Nehls, 2008), whereas other species are shown to lack the genetic capacity to act as saprotrophs (Lindahl and Tunlid, 2015; Frey, 2019).

Because of the knowledge gaps regarding the actual role of mycorrhizal fungi functioning in ecosystems, they are rarely represented in state-of-the-art TBMs that are used to simulate the global carbon and nutrient cycle responses to global change (Le Quéré et al., 2018). Most TBMs predict a higher N limitation effect on CO₂-fertilised growth than observed, potentially due to a significant underestimation of plant N acquisition enhancement under eCO₂ (Zaehle et al., 2014). One possible reason for this is that most TBMs represent N acquisition as plant root uptake of mineral N only but lack a representation of a plant-controlled pathway to increase N availability and acquisition on demand (Knops et al., 2002; Chapman et al., 2006; Zaehle and Dalmonech, 2011). This could be any pathway, which may involve direct organic nutrient uptake and/or rhizodeposition in order accelerate SOM decomposition, but in particular mycorrhiza-mediated rhizosphere–SOM interactions (Chapman et al., 2006; Phillips et al., 2013; Brzostek et al., 2017). The absence of such pathways can lead to incorrect predictions of plant and ecosystem responses to eCO₂ and to uncertain estimates of the ecosystem impact on climate change (Brzostek et al., 2017; Terrer et al., 2018; Shi et al., 2019).

Consequently, several approaches to implement mycorrhizal dynamics into ecosystem models have been developed during the last decade, either to explore C–nutrient exchange/competition between the host plant and mycorrhizal fungi for individuals (e.g. Meyer et al., 2010; Franklin et al., 2014) or to improve the representation of plant nutrient acquisition in TBMs at the ecosystem scale by simulating their effects on SOM dynamics (e.g. Deckmyn et al., 2011; Orwin et al., 2011) or by including them directly into plant N nutrition calculations (e.g. Phillips et al., 2013; Brzostek et al., 2014, 2015; Sulman et al., 2017). One of the main issues facing such a model is that both species of arbuscular and ectomycorrhizal fungi have been shown to actively mine SOM as a nutrient source, either by a direct enzymatic breakdown of material, which potentially includes saprotrophic behaviour (Hodge et al., 2001; Wu et al., 2005; Pérez-Tienda et al., 2012; Hodge and Storer, 2015), or by passively priming SOM by C exudation, which accelerates SOM decomposition (Hodge et al., 2001; Cheng et al., 2012; Lindahl and Tunlid, 2015; Paterson et al., 2016; Lang et al., 2021). Additionally, associations with either AMs or ECMs are identified at the species level (and it is difficult to aggregate them to the ecosystem level) or appointed to plant functional types (PFTs) (which are represented in TBMs). Thus, for the purpose of this study, we choose not to implement the common mycorrhizal types; rather, we explore what the consequences of these two mycorrhizal functionalities (i.e. an active SOM-mining saprotrophic mycorrhizal function and a passive SOM-priming mycorrhizal function that accelerate SOM decomposition by C exudation) are on plant N acquisition and ecosystem carbon storage, how explicitly we can model plant–mycorrhiza–soil interactions given the current knowledge and data availability, and, finally, whether or not this improves the prediction of the observed responses of plant N acquisition at the ecosystem scale under ambient- and elevated-CO₂ conditions.

We implement the two alternative representations of mycorrhizal functionalities into the QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system; Thum et al., 2019) model and then compare the model predictions for these representations to the standard model, which does not account for mycorrhizal fungi, to assess the robustness of the alternative model structures and their parameterisation. We test all model variants at the Duke Free-air CO₂ Enrichment (FACE) site to especially investigate responses to atmospheric CO₂ enhancement and, subsequently, compare model output to observational data, where available (McCarthy et al., 2006).

2 Methods

2.1 QUINCY model

The QUINCY TBM is described in detail in Thum et al. (2019). The model simulates photosynthesis and the growth of vegetation represented by average individuals of PFTs as well as soil biogeochemical processes. QUINCY represents fully coupled C, N, P and water cycles, and it traces the isotopic signals of ¹³C, ¹⁴C and ¹⁵N. QUINCY operates on a half-hourly time step to account for short-term acclimation, but many processes have a specific memory time, e.g. days, weeks or months, over which the calculated flux is averaged to also consider longer response times. Below, we briefly present the representation of the most relevant processes for plant N acquisition as well as of soil organic matter dynamics that are affected or modified by the implementation of the mycorrhizal processes.

Plant N acquisition is described as direct uptake of inorganic nitrogen by fine roots. Potential root uptake is assumed to be a function of fine-root biomass and mineral soil N content, i.e. ammonium (NH₄) and nitrate (NO₃). This potential uptake is scaled down in times of low plant N demand, which depends on the ratio of available to required N for growth. Biomass is allocated dynamically above- and below-ground in response to plant nutrient and water demand. Acquired N and assimilated C are first added to the plant labile pool, which is one of two nonstructural plant pools represented by QUINCY. From there, C and N are distributed among all other plant pools, i.e. the plant structural pool to represent tissue growth and the plant reserve pool as second nonstructural pool for storage, within a week. This allows QUINCY to account for sink and source limitation and also smooths daily variations.

QUINCY includes a representation of soil carbon and nutrient pools that is based on the CENTURY soil model (Parton et al., 1993). It accounts for five organic soil pools (i.e. three litter pools, a fast-overturning SOM pool and a slow-overturning SOM pool) and three inorganic soil pools (i.e. NH₄, NO₃ and PO₄). Decomposition is described as first-order decay by specific turnover times and pool sizes. Net mineralisation of nutrients occurs if the nutrient immobilisation for litter decomposition is less than the gross nutrient mineralisation from the decomposition of SOM. Plant uptake competes with immobilisation as well as leaching and nitrification- and denitrification-related gaseous losses for the available inorganic nutrients.

2.2 MYC model description

The plant–MYCorrhiza–soil interaction (MYC) model includes plant–mycorrhiza dynamics and mycorrhiza–soil interactions. Relevant QUINCY fluxes and MYC model fluxes are presented in Fig. 1.

Figure 1 Schematic representation of the plant–MYCorrhiza–soil interaction (MYC) model fluxes within the QUINCY framework. Boxes represent model pools: light green represents C pools and red represents N pools. They are grouped for plants (dark green), mycorrhizae (yellow) and SOM (brown). Arrows indicate fluxes: solid arrows are explained in Eqs. (1)–(12), whereas dashed arrows are not explained in detail here, as the latter are not modified by or used within the MYC model. Grey arrows belong to the QUINCY framework, black arrows belong to both mycorrhizal functionalities, blue arrows belong only to saprotrophic mycorrhizae, and orange arrows belong to or are affected by decomposing mycorrhizae. E depicts exchange fluxes, R denotes respiration fluxes, T represents turnover fluxes and U depicts uptake fluxes. Subscripts indicate whether plants (p), mycorrhizae (m) or SOM determine flux rates. Subscripts contain additional information, e.g. exchanged property or source. For more details, see Eqs. (1)–(12).

Plants and mycorrhizae interact through the C flux from host plants to mycorrhizal fungi, the reduction in the fine-root uptake capacity in the presence of mycorrhizae (as these roots are assumed to be partially covered by the mycelium and therefore have a less direct contact with the soil), and the export of nitrogen taken up by mycorrhizae to plant roots. The respective net C and N acquisition of plants (Fig. 2a, b) is described as follows:

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\mathrm{d}{C}^{\mathrm{p}}}{\mathrm{d}t}}=\mathrm{BP}=\mathrm{NPP}+\mathrm{\Delta }{E}_{\mathrm{C}},\text{(1b)}& {\displaystyle \frac{\mathrm{d}{N}^{\mathrm{p}}}{\mathrm{d}t}}=A_{\mathrm{N}}^{\mathrm{p}}=U_{\mathrm{N}}^{\mathrm{p}}+\mathrm{\Delta }{E}_{\mathrm{N}}.\end{array}$$

The net C acquisition of plants, i.e. plant biomass production (BP), is given by net primary production (NPP) and net C exchange between the host plant and mycorrhizal fungi (ΔE_C in Fig. 2a). Plant (net) N acquisition ($A_{\mathrm{N}}^{\mathrm{p}}$) is given by direct root uptake ($U_{\mathrm{N}}^{\mathrm{p}}$) and net N exchange (ΔE_N in Fig. 2b).

NPP and root N uptake are calculated by the following expressions in the standard QUINCY model, which have been further modified to account for the effects of mycorrhizae:

$$\begin{array}{}\text{(2a)}& {\displaystyle }\begin{array}{rl}{U}_{\mathrm{N},j}^{\mathrm{p},\mathrm{pot},*}& =v_{max,j}^{*}\times N_j\times \left(K_{\mathrm{m}\mathrm{1},j}+{\displaystyle \frac{\mathrm{1}}{K_{\mathrm{m}\mathrm{2},j}+N_j}}\right)\\ & \times {\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}}\times C_{\mathrm{fr}}^{\mathrm{p},*},\text{where}\end{array}\text{(2b)}& {\displaystyle }{v}_{max,j}^{*}={\displaystyle \frac{v_{max,j}}{\mathrm{1}-f_{\mathrm{cover},\max }+f_{\mathrm{cover},\max }\times f_{U_{\mathrm{N}}^{\mathrm{eff}}}^{\mathrm{m}}}}\text{and}\text{(2c)}& {\displaystyle }{C}_{\mathrm{fr}}^{\mathrm{p},*}=C_{\mathrm{fr}}^{\mathrm{p}}\times f_{\mathrm{cover}}.\text{Here,}\text{(2d)}& {\displaystyle }{f}_{\mathrm{cover}}=\mathrm{MIN}\left(f_{\mathrm{cover},\max },{\displaystyle \frac{C^{\mathrm{m}}}{C_{\mathrm{fr}}^{\mathrm{p}}}}\right).\end{array}$$

In the above expressions, j refers to either NH₄ or NO₃, N_j denotes the NH₄ or NO₃ concentration in soil, $v_{max,N_j}^{*}$ is a parameter describing the maximum uptake rate per unit fine-root biomass ($C_{\mathrm{fr}}^{\mathrm{p},*}$), $K_{\mathrm{m}\mathrm{1},N_j}$ and $K_{\mathrm{m}\mathrm{2},N_j}$ are the respective low- and high-affinity half-saturation parameters, and ${\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}}$ is a scaling factor that describes current plant N demand. The modification of the maximum uptake rate per unit biomass ($v_{max,j}^{*}$) maintains the overall maximum uptake capacity of the rhizosphere ($v_{max,j}$) by taking the relative presence of mycorrhizae (f_cover,max) and their higher uptake efficiency ($f_{U_{\mathrm{N}}^{\mathrm{eff}}}^{\mathrm{m}}$) compared to plants into account, while the adjustment of root biomass considers the disconnection of root tips from soil by mycorrhizal fungi. This coverage is given by the current ratio of fine-root biomass ($C_{\mathrm{fr}}^{\mathrm{p}}$) and mycorrhizal biomass (C^m), but it is constrained by a maximum coverage to account for a shift in the mycorrhizal community to more exploitative species if plants are under severe N limitation.

In a similar way to plant C and N acquisition (Eq. 1), mycorrhizal growth (G^m in Fig. 2c, d) is defined as the net budget for C and N accumulation over time. $G_{\mathrm{C}}^{\mathrm{m}}$ is given by mycorrhizal C uptake from SOM ($U_{\mathrm{C}}^{\mathrm{m}}$, only in the case of saprotrophic mycorrhizae); net C exchange between the host plant and mycorrhizal fungi (ΔE_C); and losses by mycorrhizal turnover ($T_{\mathrm{C}}^{\mathrm{m}}$), respiration (R^m) and exudation to accelerate SOM decomposition ($D_{\mathrm{SOM}}^{\mathrm{m}}$, only in the case of decomposing mycorrhizae; Fig. 2c). $G_{\mathrm{N}}^{\mathrm{m}}$ is given by mycorrhizal N uptake ($U_{\mathrm{N}}^{\mathrm{m}}$), net N exchange (ΔE_N) and mycorrhizal turnover ($T_{\mathrm{N}}^{\mathrm{m}}$ in Fig. 2d).

$$\begin{array}{}\text{(3a)}& {\displaystyle \frac{\mathrm{d}{C}^{\mathrm{m}}}{\mathrm{d}t}}=G_{\mathrm{C}}^{\mathrm{m}}=U_{\mathrm{C}}^{\mathrm{m}}-\mathrm{\Delta }{E}_{\mathrm{C}}-T_{\mathrm{C}}^{\mathrm{m}}-R^{\mathrm{m}}-D_{\mathrm{SOM}}^{\mathrm{m}}\text{(3b)}& {\displaystyle \frac{\mathrm{d}{N}^{\mathrm{m}}}{\mathrm{d}t}}=G_{\mathrm{N}}^{\mathrm{m}}=U_{\mathrm{N}}^{\mathrm{m}}-\mathrm{\Delta }{E}_{\mathrm{N}}-T_{\mathrm{N}}^{\mathrm{m}}\end{array}$$

Figure 2 C and N budgets of plants (a, b) and mycorrhizal fungi (c, d), as described in Eqs. (1) and (3), per day (of year, DOY). Net budgets are shown in red, C fluxes are given in black and N fluxes are presented in grey. Uptake (U) is shown using dashed lines, net exchange rates (ΔE) are presented using solid lines, losses by turnover (T) and respiration (R) are shown using dot-dash lines (c, d), and losses by exudation to SOM (D) are presented using long dashes (c).

The exchange of C and N between the host plant and mycorrhizal fungi consists of an exudation flux from plants to mycorrhizae (E^p2m), where we assume that plant N exudation is zero and that the export flux from mycorrhizae to plants (E^m2p in Fig. 3) is as follows:

$$\begin{array}{}\text{(4a)}& {\displaystyle }\mathrm{\Delta }{E}_{\mathrm{C}}=+E_{\mathrm{C}}^{\mathrm{m}\mathrm{2}\mathrm{p}}-E_{\mathrm{C}}^{\mathrm{p}\mathrm{2}\mathrm{m}},\text{(4b)}& {\displaystyle }\mathrm{\Delta }{E}_{\mathrm{N}}=+E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}.\end{array}$$

Positive signs depict fluxes from mycorrhizae to plants, whereas negative signs depict fluxes from plants to mycorrhizae.

Figure 3 C and N exchange rates between host plants and mycorrhizal fungi, as described in Eq. (4), per day (of year, DOY). C fluxes are shown in black and N fluxes are given in grey, as in Fig. 2. Net exchange rates (ΔE, solid lines) are partitioned into exudation from plants to mycorrhizal fungi (E^p2m, dashed line) and export from mycorrhizal fungi to plants (E^m2p, dot-dash line).

We constrain C exudation ($E_{\mathrm{C}}^{\mathrm{p}\mathrm{2}\mathrm{m}}$) from the host plant to mycorrhizal fungi in a similar way to that assumed in MYCOFON using the following parameters (Meyer et al., 2010):

• the amount of available C, which is represented by the labile C pool in QUINCY;

• plant N demand, where increasing N demand increases C allocation to mycorrhizal growth;

• an interaction term that allows for a decline in exudation if plants do not receive nutrients in return, which could happen either in the case of satisfied demand by plants or in the case of severe N limitation, thereby turning mycorrhizal fungi into competitors that do not deliver any N to host plants;

• a minimum and maximum amount of mycorrhizae compared with plant fine roots, to avoid complete mycorrhizal death and to prevent plants overspending C.

For a detailed description of the equations, see Appendix A.

We calculate N export from mycorrhizal fungi to host plants ($E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}$) as mycorrhizal surplus N, which is derived from current mycorrhizal N uptake ($U_{\mathrm{N}}^{\mathrm{m}}$) and N demand (${\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{m}}$), assuming that mycorrhizal fungi will first maintain their own C : N ratio (${\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}}$) before they can deliver N to host plants. This potential N export is scaled by plant N demand (${\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}}$) for actual N export to suppress export if the plants have sufficient N. As we assume that mycorrhizal fungi deliver N as part of amino acids, C export ($E_{\mathrm{C}}^{\mathrm{m}\mathrm{2}\mathrm{p}}$) is given by N export and a prescribed C : N ratio (${\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}\mathrm{2}\mathrm{p}}$).

$$\begin{array}{}\text{(5a)}& {\displaystyle }{E}_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}=(U_{\mathrm{N}}^{\mathrm{m}}-{\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{m}})\times {\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}},\text{with}\text{(5b)}& {\displaystyle }{\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{m}}={\displaystyle \frac{C^{\mathrm{m}}}{{\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}}}}-N^{\mathrm{m}},\text{(5c)}& {\displaystyle }{E}_{\mathrm{C}}^{\mathrm{m}\mathrm{2}\mathrm{p}}={\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}\mathrm{2}\mathrm{p}}\times E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}.\end{array}$$

Mycorrhizal N uptake ($U_{\mathrm{N}}^{\mathrm{m}}$) is the sum of mineral N uptake ($U_{N_j}^{\mathrm{m}}$) and organic N uptake ($U_{N_{\mathrm{SOM}}}^{\mathrm{m}}$) in the case of saprotrophic mycorrhizae. Mineral N uptake, i.e. uptake of ammonium (NH₄) and nitrate (NO₃), is calculated similarly to plant root uptake in QUINCY (Eq. 2), but it is down-regulated if the C investment in N transformation to amino acids exceeds a threshold to avoid mycorrhizal death:

$$\begin{array}{}\text{(6a)}& {\displaystyle }\begin{array}{rl}{U}_{N_j}^{\mathrm{m},\mathrm{pot}}& =v_{max,N_j}^{\mathrm{m}}\times N_j\times \left(K_{\mathrm{m}\mathrm{1},N_j}+{\displaystyle \frac{\mathrm{1}}{K_{\mathrm{m}\mathrm{2},N_j}+N_j}}\right)\\ & \times C^{\mathrm{m}}\times f_{\mathrm{invest}},\text{where}\end{array}\text{(6b)}& {\displaystyle }{f}_{\mathrm{invest}}={\displaystyle \frac{f_{\mathrm{resp},\max }\times C^{\mathrm{m}}}{R_{U_{N_j}}^{\mathrm{m},\mathrm{pot}}}}.\end{array}$$

Here, N_j denotes either NH₄ or NO₃, $v_{max,N_j}^{\mathrm{m}}$ is a parameter describing the maximum uptake rate per unit mycorrhizal biomass (C^m), $K_{\mathrm{m}\mathrm{1},N_j}$ and $K_{\mathrm{m}\mathrm{2},N_j}$ are the respective low- and high-affinity half-saturation parameters, and f_resp,max is the empirical threshold for C investment into mineral N uptake.

We calculate the turnover of mycorrhizal fungi ($T_X^{\mathrm{m}}$) based on the life span of mycorrhizae as follows:

$$\begin{array}{}\text{(7)}& T_X^{\mathrm{m}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{myc}}}}\times X^{\mathrm{m}},\end{array}$$

where τ_myc is mycorrhizal turnover time and X^m is either mycorrhizal C or N. $T_X^{\mathrm{m}}$ is added to the respective litter pool in soil, similar to the death plant fine roots.

We calculate respiration of mycorrhizal fungi (R^m) as the sum of direct mycorrhizal respiration ($R_{\mathrm{gm}}^{\mathrm{m}}$), which is associated with growth and maintenance processes, using a constant C use efficiency (CUE^m) and recently gained C, and N uptake respiration ($R_{U_{N_j}}^{\mathrm{m}}$), which is related to the specific costs (r_j) associated with transforming mineral N (N_j) into amino acids (Zerihun et al., 1998).

$$\begin{array}{}\text{(8a)}& {\displaystyle }{R}^{\mathrm{m}}=R_{\mathrm{gm}}^{\mathrm{m}}+R_{U_{N_j}}^{\mathrm{m}}\text{(8b)}& {\displaystyle }{R}_{\mathrm{gm}}^{\mathrm{m}}=(\mathrm{1}-{\mathrm{CUE}}^{\mathrm{m}})\times (U_{\mathrm{C}}^{\mathrm{m}}+E_{\mathrm{C}}^{\mathrm{p}\mathrm{2}\mathrm{m}})\text{(8c)}& {\displaystyle }{R}_{U_{N_j}}^{\mathrm{m}}=r_j\times U_{N_j}^{\mathrm{m}}\end{array}$$

Special functionalities of saprotrophic and decomposing mycorrhizae

Saprotrophic mycorrhizae mine SOM as additional organic N source. We assume that mycorrhizae access old, C-depleted, N-rich SOM. We further assume that mycorrhizae are only facultative saprotrophs, which is why we decrease N uptake from SOM if both plants and mycorrhizae have sufficient N, which is quantified by the current C : N ratio of mycorrhizae.

$$\begin{array}{}\text{(9a)}& {\displaystyle }\begin{array}{rl}& U_{N_{\mathrm{SOM}}}^{\mathrm{m}}=v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}{\displaystyle \frac{N_{\mathrm{SOM}}}{K_{\mathrm{m},N_{\mathrm{SOM}}}+N_{\mathrm{SOM}}}}\times C^{\mathrm{m}}\times f_{\mathrm{limit}},\\ & \text{where}\end{array}\text{(9b)}& {\displaystyle }{f}_{\mathrm{limit}}=\left\{\begin{array}{ll}\mathrm{1},& \text{if}\frac{C^{\mathrm{m}}}{N^{\mathrm{m}}}>{\mathit{\chi }}_{\mathrm{CN},\frac{min}{\mathrm{2}}}^{\mathrm{m}}\\ \mathrm{0.5},& \text{if}{\mathit{\chi }}_{\mathrm{CN},\frac{min}{\mathrm{2}}}^{\mathrm{m}}>\frac{C^{\mathrm{m}}}{N^{\mathrm{m}}}>{\mathit{\chi }}_{\mathrm{CN},min}^{\mathrm{m}},\\ & \text{with}{\mathit{\chi }}_{\mathrm{CN},\frac{min}{\mathrm{2}}}^{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{2}}\times ({\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}}+{\mathit{\chi }}_{\mathrm{CN},min}^{\mathrm{m}})\\ \mathrm{0},& \text{if}\frac{C^{\mathrm{m}}}{N^{\mathrm{m}}}<{\mathit{\chi }}_{\mathrm{CN},min}^{\mathrm{m}}.\end{array}\right.\end{array}$$

Here, $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$ is the maximum uptake rate per unit mycorrhizal biomass (C^m); $K_{\mathrm{m},N_{\mathrm{SOM}}}$ is the half-saturation parameter; and f_limit is a down-regulation function for N uptake from SOM, which is derived from the target C : N ratio of mycorrhizal fungi (${\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{m}}$), a minimum C : N ratio (${\mathit{\chi }}_{\mathrm{CN},\min }^{\mathrm{m}}$) and the current C : N ratio.

We assume that organic N uptake by SOM-mining mycorrhizae is associated with a C flux corresponding to the uptake of amino acids and some additional C, even if most mycorrhizal fungi lack the capacity to assimilate soil organic C (SOC) and acquire C only from their host plant (Lindahl and Tunlid, 2015; Frey, 2019). This C uptake ($U_{\mathrm{C}}^{\mathrm{m}}$) is described as follows:

$$\begin{array}{}\text{(10)}& U_{\mathrm{C}}^{\mathrm{m}}=U_{N_{\mathrm{SOM}}}^{\mathrm{m}}\times {\mathit{\chi }}_{\mathrm{CN}}^{\mathrm{SOM}}\times {\mathrm{CUE}}^{\mathrm{m}},\end{array}$$

where CUE^m is the C use efficiency of mycorrhizal uptake and describes the share of SOC that mycorrhizae can incorporate. The remaining C is respired.

Decomposing mycorrhizae exude C ($D_{\mathrm{SOM}}^{\mathrm{m}}$) to accelerate decomposition. This C provides energy for microbes to decompose older SOM faster, which enhances the net mineralisation flux and, therefore, increases soil mineral N availability for mycorrhizal fungi as well as plants and microbes. We calculate the exudation as follows:

$$\begin{array}{}\text{(11)}& D_{\mathrm{SOM}}^{\mathrm{m}}=f_{D_{\mathrm{SOM}}}\times C^{\mathrm{m}},\end{array}$$

where $f_{D_{\mathrm{SOM}}}$ is an empirical parameter. In order to not violate the fixed soil C : N ratio assumed in QUINCY, $D_{\mathrm{SOM}}^{\mathrm{m}}$ is added to heterotrophic respiration. This is reasonable, as we assume that the induced acceleration of decomposition increases heterotrophic respiration generally.

We represent the resulting acceleration of decomposition implicitly by lowering the turnover time of slow SOM in the presence of decomposing mycorrhizal fungi.

$$\begin{array}{}\text{(12)}& {\mathit{\tau }}_{\mathrm{SOM}}^{*+}={\mathit{\tau }}_{\mathrm{SOM}}^{*}\times {\displaystyle \frac{K_{\mathrm{m},{\mathit{\tau }}_{\mathrm{SOM}}}+C^{\mathrm{m}}}{v_{max,{\mathit{\tau }}_{\mathrm{SOM}}}\times C^{\mathrm{m}}}},\end{array}$$

where ${\mathit{\tau }}_{\mathrm{SOM}}^{*}$ denotes the turnover time of the SOM without the influence of decomposers but in mycorrhizal presence, $K_{\mathrm{m},{\mathit{\tau }}_{\mathrm{SOM}}}$ and $v_{max,{\mathit{\tau }}_{\mathrm{SOM}}}$ are empirical parameters, and C^m is mycorrhizal biomass.

In order to maintain the overall ecosystem dynamics in a steady state, comparable to the model simulations without mycorrhizae, we multiply the turnover of the slow-overturning SOM pool in QUINCY by an empirically determined parameter $f_{{\mathit{\tau }}_{\mathrm{SOM}}}$ for both mycorrhizal functionalities.

Meyer et al. (2010)Meyer et al. (2010)Marschner and Dell (1994)Göransson et al. (2006)Marschner and Dell (1994)Göransson et al. (2006)Marschner and Dell (1994)Smith and Read (2010)Allen et al. (2003)Wallander et al. (2003)Hauptmann (1985)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)Thum et al. (2019)

Table 1 MYC model parameters.

2.3 Duke Forest FACE experiment

We test the QUINCY–MYC model at the Duke Forest Free-air CO₂ Enrichment (FACE) site (hereafter Duke FACE). The Duke FACE experiment (McCarthy et al., 2007) was set up in a loblolly pine (Pinus taeda) plantation (35.9° N, 79.08° W; North Carolina, US) that was established in 1983 after a clear-cut. The experiment started in August 1996 with three plots per CO₂ treatment, i.e. ambient and elevated (+200 ppm) CO₂ concentration, with plots paired according to soil N availability.

For our analysis, we use locally measured climate data (see Sect. 2.4), comprising annual measurement data for ambient and elevated treatments for plant C and N fluxes, i.e. plant biomass production (BP) and plant N acquisition ($A_{\mathrm{N}}^{\mathrm{p}}$), which are estimated by Finzi et al. (2007). Responses were calculated per each treatment plot pair and then averaged. Soil C concentrations and δ¹³C measurements are taken from Lichter et al. (2008), who used specific δ¹³C labelling of the artificial CO₂ that was used for fumigation.

2.4 Model set-up and protocol

QUINCY requires half-hourly meteorological forcing consisting of air temperature, precipitation (rain and snow), long-wave and short-wave radiation, atmospheric CO₂ concentration, and N and P deposition rates. Meteorological data are derived from version 7 of the CRUNCEP dataset (Viovy, 2016), which provides daily data from 1901 to 2015. Data are disaggregated using a statistical weather generator (Zaehle and Friend, 2010) to the half-hourly model time step. Annual data for atmospheric CO₂ concentration are taken from Le Quéré et al. (2018), while N deposition rates are from Lamarque et al. (2010) and Lamarque et al. (2011). As local records of forcing data are available, we used data from Walker et al. (2014) for the duration of the Duke FACE experiment.

Additionally, QUINCY requires information about vegetation, given as PFT and soil information, such as texture, bulk density and rooting and soil depth, for each site. As this work aims to study N controls on plant growth, we also use information about inorganic P content that is kept constant to avoid the influence of P limitation.

QUINCY pools (vegetation and soil) are brought to quasi-equilibrium by a spin-up period of 500 years that uses repeated meteorological data from 1901 to 1930 to drive simulations, before the actual simulation period (1901–2015) starts, with transient climate and CO₂ concentrations as described above. Vegetation tissue pools are set to zero in 1983, the year of forest establishment at the site. Harvested biomass persists in the system by adding it to litter, except for woody biomass, from which a fraction of 80 % is removed.

We perform two simulations: one for ambient CO₂ and one in which atmospheric CO₂ concentrations are elevated by 200 ppm from August 1996 on in order to mirror the FACE experiment. The δ¹³C signal is calculated by the amount of fumigated artificial CO₂ compared with natural CO₂ (McCarthy et al., 2007).

To assess parameter uncertainty, we vary all MYC model parameters (Table 1, section above the line) within a range of ±10 %. We use Latin hypercube sampling (LHS; Saltelli et al., 2000) to generate 200 parameter sets for each mycorrhizal functionality. We did not vary default QUINCY parameters (Table 1, section below the line), as a detailed analysis can be found in Thum et al. (2019).

For the analysis of the Duke FACE experiment, we carry out two runs for each model variant and each parameter set: QUINCY without mycorrhizal fungi (1 times two runs, i.e. one under ambient CO₂ and one under elevated CO₂, as described above), QUINCY with decomposing mycorrhizae (200 times two runs, i.e. one run for each CO₂ concentration for each of the 200 LHS parameter sets) and QUINCY with saprotrophic mycorrhizae (200 times two runs, i.e. one run for each CO₂ concentration for each of the 200 LHS parameter sets). We analyse the median and the 10 % and 90 % quantiles of the LHS simulations.

In order to better present changes in the CO₂ elevation during the experiment, flux responses are shown as relative responses (δY in Eq. 13a), whereas changes in pool sizes, e.g. soil organic matter content or mycorrhizal biomass, are shown as absolute changes (ΔX in Eq. 13b).

$$\begin{array}{}\text{(13a)}& {\displaystyle }\mathit{\delta }Y={\displaystyle \frac{Y_{\mathrm{e}}-Y_{\mathrm{a}}}{Y_{\mathrm{a}}}}\text{(13b)}& {\displaystyle }\mathrm{\Delta }X=X_{\mathrm{e}}-X_{\mathrm{a}}\end{array}$$

In order to better understand the impact of the mycorrhizal parameterisation on the simulated ecosystem dynamics, we also use a linear decomposition method (Rastetter et al., 1992) to attribute the changes in ecosystem C storage under eCO₂ to changes in the N cycle induced by mycorrhizae.

3 Results

3.1 Effects on simulated vegetation dynamics under ambient and elevated CO₂

Simulated ambient plant biomass production (BP) is generally lower than observed BP (999.6 ± 138.8 gC m⁻² yr⁻¹), but it does not differ much among simulations with and without mycorrhizal fungi (Fig. 4a). Simulated ambient plant N acquisition ($A_{\mathrm{N}}^{\mathrm{p}}$) captures observed plant N acquisition (8.3 ± 0.5 gN m⁻² yr⁻¹) in general (Fig. 4c), whereby both BP and N acquisition are higher in simulations without mycorrhizae (BP: 808.7 ± 39.5 gC m⁻² yr⁻¹; $A_{\mathrm{N}}^{\mathrm{p}}$: 9.0 ± 0.2 gN m⁻² yr⁻¹) and with saprotrophic mycorrhizae (BP: 755.6 ± 41.1 gC m⁻² yr⁻¹; $A_{\mathrm{N}}^{\mathrm{p}}$: 9.8 ± 1.4 gN m⁻² yr⁻¹), whereas these values are lower in simulations with decomposing mycorrhizae (BP: 678.1 ± 71.2 gC m⁻² yr⁻¹; $A_{\mathrm{N}}^{\mathrm{p}}$: 6.3 ± 0.5 gN m⁻² yr⁻¹). This indicates a stronger competition for N by decomposing mycorrhizae than by saprotrophic mycorrhizae, which in turn limits plant growth. Root infection by mycorrhizal fungi lowers direct plant N uptake, as mycorrhizae reduce root contact with the soil ($U_{\mathrm{N}}^{\mathrm{p}}$ in Fig. 4e) by around 84 %; however, this effect is balanced the by mycorrhizal N uptake and export to plants (ΔE_N in Fig. 4g). N export is higher for saprotrophic mycorrhizae (8.4 ± 1.4 gN m⁻² yr⁻¹) than for decomposing mycorrhizae (4.5 ± 0.4 gN m⁻² yr⁻¹), leading to the simulated difference in total plant N acquisition.

Elevated CO₂ (eCO₂) increases both observed and modelled BP. The first-year effect is lower for observations (+24.8 %) and for simulations with saprotrophic mycorrhizae (+25.0 %), whereas it is higher for simulations without mycorrhizae (+43.6 %) and with decomposing mycorrhizae (+35.3 %; δBP in Fig. 4b). The observed BP response reaches a value of +27.2 ± 2.1 % for the first 6 years of the experiment; it then increases to +36.5 ± 2.9 % for 5 years and decreases slightly in the last year (+20.3 %). This behaviour is not entirely captured by any of the simulations, but it also remains unclear whether the last-year decrease is a sign of progressive N limitation caused by annual variation or by a bias in observations after winter 2001–2002, when a storm affected some of the plots (McCarthy et al., 2006). The positive observed responses of plant N acquisition during the entire experimental phase (+23.5 ± 6.2 %; $\mathit{\delta }{A}_{\mathrm{N}}^{\mathrm{p}}$ in Fig. 4c) point to the effect of climate variability and disturbance, rather than progressive N limitation (McCarthy et al., 2006). Modelled responses of BP are similar in simulations without mycorrhizal fungi and with saprotrophic mycorrhizae within the first half of the experiment. In the second half of the experiment, the BP response simulated by the model without mycorrhizal fungi decreases. This is constrained by the small increase in simulated plant N acquisition within the first 5 years (+8.9 ± 4.0 %) and the slightly negative response afterwards (−5.0 ± 5.5 %). In contrast, the simulated BP response of the model with saprotrophic mycorrhizae remains enhanced until the end of the experiment (+39.6 ± 1.2 %), which is supported by a strong positive response of plant N acquisition during the entire experiment (+33.8 ± 12.9 %). The simulated BP response by the model variant with decomposing mycorrhizae is only positive during the first 3 years of the experiment and is slightly negative on average (−3.3 ± 16.1 %). This response is caused by a small but entirely negative response of plant N acquisition (−2.6 ± 1.2 %). Notably, this negative response of plant N acquisition is not caused by reduced plant N uptake in response to eCO₂ as in the model variant without mycorrhizal fungi ($\mathit{\delta }{U}_{\mathrm{N}}^{\mathrm{p}}$ in Fig. 4f). In the simulations, decomposing mycorrhizae accelerate N decomposition in response to eCO₂ and, therefore, allow plants to increase N uptake (+6.8 ± 2.8 %). However, they also act as strong competitors and, while mycorrhizae take up N from the soil, they reduce N export to their host plants to satisfy their own N demand (−7.6 ± 2.8 %; δΔE_N in Fig. 4h). In contrast, the simulated positive plant N acquisition response during the entire experiment in the model variant with saprotrophic mycorrhizae is caused by both a slight increase in plant N uptake (+1.8 ± 5.7 %) and a strong increase in mycorrhizal N export (+40.3 ± 14.9 %), which matches the observed total plant N acquisition response best.

Figure 4 Vegetation C and N fluxes for the Duke FACE experiment for plant biomass production (BP; a, b), plant N acquisition ($A_{\mathrm{N}}^{\mathrm{p}}$; c, d), plant N uptake ($U_{\mathrm{N}}^{\mathrm{p}}$; e, f) and N export from mycorrhizae to plant ($E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}$; g, h) from ambient-CO₂ simulations and treatments (a, c, e, g) and in response to elevated CO₂ (b, d, f, h), which is calculated by Eq. (13a). The simulations are presented as follows: grey – QUINCY without mycorrhizae; orange – QUINCY with decomposing mycorrhizae; blue – QUINCY with saprotrophic mycorrhizae. The line represents the LHS simulation median and shaded areas present the 10 %–90 % quantile range of LHS simulations. Black denotes observations by Finzi et al. (2007) averaged over treatments, with the error bars indicating ±1 SE.

Table 2 Change in mycorrhizal growth during the Duke FACE experiment.

Differences in simulated N export rates from mycorrhizae to host plants (ΔE_N) under ambient CO₂ among mycorrhizal types are caused by differences in total N uptake by mycorrhizae ($U_{\mathrm{N}}^{\mathrm{m}}$ in Fig. 5a). Generally, decomposing mycorrhizae take up almost 30 % more N from mineral sources than saprotrophic mycorrhizae ($U_{N_{min}}^{\mathrm{m}}$ in Fig. 5c; decomposers: 8.1 ± 0.4 gN m⁻² yr⁻¹; saprotrophs: 6.3 ± 0.8 gN m⁻² yr⁻¹). However accessing SOM as an N source ($U_{N_{\mathrm{SOM}}}^{\mathrm{m}}$ in Fig. 5e; 7.2 ± 1.5 gN m⁻² yr⁻¹) allows saprotrophic mycorrhizae to more than double their N uptake, leading to a 66 % higher ambient mycorrhizal N uptake by saprotrophic mycorrhizae than by decomposing mycorrhizae.

Under eCO₂, only saprotrophic mycorrhizae can increase their mycorrhizal N uptake significantly (by +23.5 ± 8.5 %), whereas the simulated mycorrhizal N uptake response ($\mathit{\delta }{U}_{\mathrm{N}}^{\mathrm{m}}$) is almost zero for the model variant with decomposing mycorrhizae (+1.5 ± 1.5 %; Fig. 5b). This is caused by the greater competition for mineral N under eCO₂ by plants and decomposing mycorrhizal fungi, which constrains the eCO₂ response of mycorrhizal N uptake from mineral sources by decomposing mycorrhizae and plant N uptake (+4.9 ± 1.6 %) respectively (Figs. 4f, 5d). However, the slight increase in mineral N uptake in response to eCO₂, in comparison with the reduced N uptake in the model variant without mycorrhizal fungi, shows the effect of decomposing mycorrhizae on N cycling. Only the strong competition for mineral N then leads to the stronger limitation of plant growth. In contrast, saprotrophic mycorrhizae increase the uptake of both mineral (+5.5 ± 5.8 %) and organic N (+39.8 ± 13.8 %) in response to eCO₂ (Fig. 5d, f), resulting in the strong positive overall response of mycorrhizal N uptake that allows saprotrophic mycorrhizae to deliver more N to their host plants (Fig. 4h). The positive response of mycorrhizal N uptake from mineral sources also indicates an increase in mycorrhizal biomass, especially at the beginning of the experiment, which is caused by increased C exudation from plants to mycorrhizal fungi in response to eCO₂.

Figure 5 Mycorrhizal N uptake simulations for the Duke FACE experiment. Total uptake ($U_{\mathrm{N}}^{\mathrm{m}}$; a, b), uptake of mineral N ($U_{N_{\min }}^{\mathrm{m}}$; c, d), and N uptake from SOM ($U_{N_{\mathrm{SOM}}}^{\mathrm{m}}$; e, f) from ambient-CO₂ simulations (a, c, e) and in response to elevated CO₂ (b, d, f), which is calculated by Eq. (13a). The simulations are presented as follows: orange – QUINCY with decomposing mycorrhizae; blue – QUINCY with saprotrophic mycorrhizae. The line represents the LHS simulation median and shaded areas present the 10 %–90 % quantile range of LHS simulations.

Generally, the model behaviour does not change strongly for either mycorrhizal type if the MYC model parameters are varied by ±10 %. For all analysed time series and responses to eCO₂, we find a larger spread among parameter sensitivity simulations for the saprotrophic mycorrhizae than for decomposing mycorrhizae (shaded areas in Figs. 4 and 5). The N efficiency of exudation (i.e. the N gain per unit C that is exuded to mycorrhizal fungi) of saprotrophic mycorrhizae is strongly dependent on their ability to mine N from SOM (Eq. 9), which (in the model) is primarily controlled by the parameter $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$ that determines the maximum capacity of organic N uptake per unit mycorrhizal biomass. Increasing $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$ can increase the mycorrhizal surplus of N and, thus, the N export of saprotrophic mycorrhizae to plants. This creates a feedback in the model in which an increased return of N increases the amount of C that plants invest in exudation and subsequently also mycorrhizal biomass. Notably, the effect on plant C exudation can be positive (because plants support supportive mycorrhizal fungi) or negative (if mycorrhizal fungi already export enough N to fulfil plant requirements). In contrast, decomposing mycorrhizae always act as competitors for mineral N (Figs. 4e, 5c); thus, plants exude the minimum amount of C required for their N demand. This implies that the modelled magnitude of decomposing mycorrhizae and their effect on plant N nutrition is insensitive to mycorrhizal parameters but mostly controlled by the dynamically varying plant demand.

3.2 Effects on simulated soil dynamics under ambient and elevated CO₂

Mycorrhizal fungi also affect soil processes and soil organic C (SOC) stocks, either by taking up and respiring C that was attached to soil organic N in the case of saprotrophic mycorrhizae or by accelerating SOM decomposition in the case of priming mycorrhizae; this caused differences in the amount of SOC and its age during the FACE experiment. Thus, we explore changes in total SOC as well as the δ¹³C signal of SOM. This gives an indication of (1) how much of the SOC has been accumulated during the experiment in total and (2) the turnover of accumulated C, because a decrease in the δ¹³C signal of SOM indicates the respiration of younger C and, thus, potentially a faster ecosystem C turnover, as it carries the more negative δ¹³C signal of artificial fumigation CO₂ (Fig. 6a, b).

Figure 6 Soil C and N dynamics for the Duke FACE experiment. Panels (a) and (b) represent C forcing by the atmospheric CO₂ concentration (a) and related δ¹³C (b) for ambient- and elevated-CO₂ simulations. Panels (c)–(g) represent changes in soil-dynamics-related variables from the beginning of the experiment: SOC accumulation (ΔC_SOM; c) and changes in the SOM-δ¹³C (d) responses of heterotrophic respiration and the net mineralisation rate to elevated CO₂ (δR_h and δnet mineralisation in e and g respectively) and changes in R_h δ¹³C (f). Absolute changes in δ¹³C (b), SOC (c), SOM δ¹³C (d) and R_h δ¹³C (f) are calculated by Eq. (13b), whereas relative changes in heterotrophic respiration (e) and the net mineralisation rate (g) are calculated by Eq. (13a). Panels (a–d) and (f) show the ambient-CO₂ treatment (dashed line) and elevated-CO₂ treatment (solid line). Panels (c)–(g) present the simulations: grey – QUINCY without mycorrhizae; orange – QUINCY with decomposing mycorrhizae; blue – QUINCY with saprotrophic mycorrhizae. The line represents the LHS simulation median and the shaded ribbons present the 10 %–90 % quantile range of LHS simulations. Black denotes observations by Lichter et al. (2008), with the error bars indicating ±1 SE.

Observed ambient SOC (SOC_a) is increased by 187 gC m⁻² from 1996 to 2005, whereas SOC under eCO₂ (SOC_e) is decreased by 46 gC m⁻² (Fig. 6c). Additionally SOC_e has a more negative δ¹³C signal (Fig. 6d), indicating that, even though less SOC is stored, more of this SOC has been created during the experiment. We find a similar pattern for total SOC, i.e. a less strong SOC increase under eCO₂ than under aCO₂ in the model version without mycorrhizal fungi (ΔSOC_a: +136 gC m⁻² in 2005; ΔSOC_e: +38 gC m⁻²) and in the version with saprotrophic mycorrhizae (ΔSOC_a: +130 gC m⁻²; ΔSOC_e: +24 gC m⁻²; Fig. 6c). However, both models still suggest an accumulation of SOC over the experimental period, which is contrary to observations. SOC δ¹³C for the model variant without mycorrhizal fungi shows only a very small response to eCO₂ (Δδ¹³C_a: −0.1 ‰ in 2005; Δδ¹³C_e: −0.8 ‰), which indicates only minor changes in SOC turnover, whereas the δ¹³C of SOC in the model variant with saprotrophic mycorrhizae becomes significantly negative (Δδ¹³C_a: −0.5 ‰; Δδ¹³C_e: −5.6 ‰). This indicates an accelerated turnover of SOC, which is induced by SOC mining and organic N uptake (Fig. 6d). On the contrary, in the model variant with decomposing mycorrhizae, simulated total SOC evolution from 1996 to 2005 does not differ among CO₂ treatments (ΔSOC_a: +174 gC m⁻²; ΔSOC_e: +178 gC m⁻²). This is the consequence of a general tendency of this model variant to allocate below-ground C from plant exudation to mycorrhizal fungi. Although mycorrhizal fungi accelerate SOM decomposition, the effect is too small to prevent an increase in SOC_a. SOC_e development is similar to SOC_a development for the model variant with decomposing mycorrhizae, but the δ¹³C signal is lower under eCO₂ than under aCO₂ (Δδ¹³C_a: −0.3 ‰; Δδ¹³C_e: −3.7 ‰). This shows again that freshly assimilated CO₂ is transferred to the soil via exudation and mycorrhizal fungi and that the turnover of C within the ecosystem becomes faster (Fig. 6c, d).

Concurrent with the accelerated SOC turnover under eCO₂ during the experiment, heterotrophic respiration (R_h) is also increased in all model variants (Fig. 6e). In the model variant without mycorrhizal fungi, δR_h becomes positive only after 4 years, following the gradual accumulation of SOC (+2.9 ± 1.7 %). The similarity of SOC development under aCO₂ and eCO₂ gleaned by modelling decomposing mycorrhizae causes a weak response of R_h (+1.0 ± 0.8 %; Fig. 6e, g). However, the response of net N mineralisation to eCO₂ increases slightly during the experiment (2.1 ± 1.6 %), suggesting that the decomposing functionality is accelerating soil turnover. While this change is too small to have a strong effect on the vegetation, it still leads to a notable amount of additional mineral N that is taken up by both plants (+4.9 ± 1.6 %; Fig. 4f) and mycorrhizal fungi (1.5 ± 1.5 %; Fig. 5d). Simulated saprotrophic mycorrhizae access SOM as an N source and increase their organic N uptake in response to eCO₂ to satisfy their own and plant N demand (Fig. 5f). The model assumes that the mycorrhizae only take up the C associated with amino acids and then respire the remainder of the SOC in order to maintain the SOM C : N ratio (Eq. 10). This parameterisation leads to a strong increase in R_h under eCO₂ (+12.3 ± 10.1 %; Fig. 6e), associated with mycorrhizal uptake of N. The SOM decomposition unrelated to mycorrhizae is less affected, resulting in unchanged net N mineralisation (+0.5 ± 3.8 %; Fig. 6g). This is similar to the model variant without any mycorrhizae (−1.3 ± 10.6 %). However, the δ¹³C signal of R_h is significantly lower when modelling mycorrhizal fungi, resulting from the larger below-ground C flux and faster turnover of mycorrhizal biomass in the soil, in combination with the more negative δ¹³C of the freshly assimilated C (Fig. 6f).

These results are robust against parameter changes; however, as discussed in the previous section, model simulations show a larger spread among LHS simulations for saprotrophic mycorrhizae. This is due to the sensitivity of the MYC model with saprotrophic mycorrhizae to $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$, which describes the maximum C and N uptake from SOM (Sect. 3.1), and thus directly affects SOC turnover (Eq. 10).

3.3 Effects on ecosystem C storage under ambient and elevated CO₂

The linear decomposition method of Rastetter et al. (1992) allows us to better understand the impact of the mycorrhizal parameterisation on the simulated ecosystem dynamics. We attribute the changes in ecosystem C storage under eCO₂ to changes in the N cycle induced by mycorrhizae (Fig. 7), namely, (i) a change in total ecosystem N (ΔN_Eco); (ii) a change in the C : N ratio of vegetation (ΔC:N_Veg); (iii) a change in the C : N ratio of SOM (ΔC:N_SOM, including mycorrhizal fungi); and (iv) an altered partitioning of N among ecosystem compartments, i.e. between vegetation and SOM ($\mathrm{\Delta }{\mathrm{N}}_{\mathrm{Veg}}/{\mathrm{N}}_{\mathrm{Eco}}$).

We find that ΔN_Eco is only weakly changed under eCO₂ and has only minor effects on ecosystem C storage. This is because none of the model variants change ecosystem N (N_Eco) gain and the effects on N_Eco losses are minor.

For the model variant without mycorrhizal fungi, changes in ecosystem C are mainly caused by ΔC:N_Veg, which indicates a strong N limitation on plant growth, as increased C input from elevated CO₂ cannot be adequately matched by increased N uptake and the model has an insufficient capacity to increase the total below-ground C flux to reduce the excess C uptake (Fig. 7). As a result, the net C uptake response (δNPP) decreases after the initial increase (Fig. 4b).

In the presence of mycorrhizal fungi, the lower attribution of ΔC : N_Veg to a change in ecosystem C under eCO₂ indicates less N limitation on growth. In the case of the decomposing mycorrhizae, small increases in decomposition due to increased mycorrhizal growth (Table 2) result in only minor changes in soil and vegetation C : N ratios, and they are insufficient to lead to a substantial amount of N transfer from soils to vegetation, thereby leading to a small change in ecosystem C.

The improved N nutrition in the model variant with saprotrophic mycorrhizae that leads to an increase in ecosystem C is caused by a redistribution of N among ecosystem compartments, i.e. between vegetation and soil, due to the mycorrhizal mining of N from SOM and the ensuing export of this N to vegetation ($\mathrm{\Delta }{\mathrm{N}}_{\mathrm{Veg}}/{\mathrm{N}}_{\mathrm{Eco}}$ in Fig. 7). As vegetation C : N ratios are more than 3 times higher than soil C : N ratios (Kattge et al., 2011; Parton et al., 1993), a redistribution of N from soil to vegetation allows vegetation to incorporate much more additional C, thereby exceeding the relative losses that are caused by the loss of SOM N.

As for vegetation and soil compartments individually (Sect. 3.1 and 3.2), the results are robust against parameter variations; this supports the meaning of our findings, especially when focusing on ecosystem responses to eCO₂.

Figure 7 Attribution of the simulated change in ecosystem N (ΔN_Eco), vegetation C : N (ΔC : N_Veg), soil C : N (ΔC : N_SOM) and ecosystem N distribution ($\mathrm{\Delta }{\mathrm{N}}_{\mathrm{Veg}}/{\mathrm{N}}_{\mathrm{Eco}}$) to the simulated change in ecosystem C under elevated CO₂ for the model without mycorrhizal fungi (grey), with decomposing mycorrhizal fungi (orange) and with saprotrophic mycorrhizal fungi (blue). Bar heights indicate the median of the attributions and error bars indicate the 10 % and 90 % quantile spread among the LHS simulations.

4 Discussion

4.1 Representation of alternative mycorrhizal functional groups

We integrated a dynamic plant–mycorrhiza–soil (MYC) model in the QUINCY TBM to test plausible but parsimonious representations of mycorrhizal effects on ecosystem fluxes based on observational findings. To do so, we incorporated observed mycorrhizal functionalities but did not link them to mycorrhizal types, i.e. arbuscular mycorrhizae (AMs) and ectomycorrhizae (ECMs), nor to their potential ecosystem functioning, owing to contrasting observations and limited knowledge. This is contrary to other modelling studies that have differentiated between AMs and ECMs and linked their proposed effects on plant nutrition to the simulated host plants (e.g. Phillips et al., 2013; Brzostek et al., 2014; Sulman et al., 2014). The advantage of our approach is that it allows us to assess potential generic effects of mycorrhizal dynamics on simulated ecosystems independently of plant species or any prior knowledge of species associations. We also tested the effect of mycorrhizal fungi that only access mineral N in soil (results not shown), which is similar to what plant fine roots do, but the results did not differ from the model version without mycorrhizae, indicating that mycorrhiza–SOM interactions are the key to explaining their significance for plant nutrient nutrition.

We assessed the possible importance of mycorrhizal fungi to predict the future evolution of terrestrial carbon stores under elevated CO₂.

We found that only the model variant with saprotrophic mycorrhizae, in which mycorrhizae can directly access organic N bound in SOM, was able to simulate the observed growth response to eCO₂ and the underlying change in plant N acquisition consistent with observation-derived estimates from Finzi et al. (2007). The model without mycorrhizae was capable of simulating the observed increase in growth. However, this occurred only because of an overly strong increase in the C : N stoichiometry (Caldararu et al., 2020) and, thus, the nutrient use efficiency of growth, compared with simulations using other TBMs reported in Zaehle et al. (2014). This suggests that the direct access to SOM as an N source is the key feature missing in TBMs when simulating Duke FACE. This result does not provide evidence of a general mycorrhizal saprotrophic or, at least, SOM-mining behaviour, but it does demonstrate that different N uptake strategies need to be explored to improve future projections of ecosystem models. Our study indicates a strong need for both further modelling and experimental research on mycorrhizae, especially as our findings are clearly limited by the use of only one testing site. However, CO₂ experiments with this wealth of available data are rare, which is why we decided to use only the Duke FACE site and to extend our analysis to multiple observation streams available at that site, i.e. above-ground biomass observations as well as soil C measurements; due to the multiple ongoing FACE experiments in a variety of ecosystems, our analysis could be extended in the future (e.g. EucFACE, Drake et al., 2016, or BIFoR FACE, MacKenzie et al., 2021).

The model variant in which mycorrhizae do not have direct access to organic N but induce an acceleration of SOM decomposition/priming was unable to reproduce either the observed growth or N acquisition response and predicted progressive N limitation. This result occurred because, on the simulated timescale, the implicit treatment of mycorrhizal effects on SOM decomposition did not provide enough N for the plants, as mycorrhizal fungi, which are more efficient at taking up mineral N, hindered increased N export to plants because they required the additional available N to support their own N demand. This turns the decomposing mycorrhizal fungi into competitors, not facilitators, with respect to N in the short-term and, consequently, results in increased N limitation to plant growth, despite the simulated acceleration of SOM decomposition and net N mineralisation under eCO₂. However, the set-up of FACE studies, especially their rapid CO₂ enhancement, is artificial, and the observation period is rather short. Therefore, the results presented here may not represent the current environmental change as well as a full adjustment of the observed ecosystem to the change. On longer timescales, the decomposing functionality may work to improve N nutrition of plants successfully.

Furthermore, our analysis clearly showed that the exudation of freshly assimilated C to any mycorrhizal fungi is necessary to explain the observed changes in the δ¹³C signal in SOC (Lichter et al., 2008). This shows that changes in soil C stocks can only be simulated correctly when models include a form of C exudation to the soil.

4.2 MYC model simplifications and assumptions

4.2.1 MYC model limitations due to observational constraints

The MYC model contains 17 parameters. In total, 13 of these belong to both functionalities, and they control the plant–mycorrhiza interactions (e.g. C export to mycorrhiza) and general mycorrhizal dynamics (e.g. mycorrhizal stoichiometry and biomass turnover). Two parameters each describe the mycorrhizal effects on SOM specific to each functionality. We found observation-based estimates for only five parameters (Table 1). This lack of suitable data is caused by the large difficulties involved with measuring within the rhizosphere without destroying plant–mycorrhiza symbioses (Savage et al., 2018; Vicca et al., 2018; Lang et al., 2021), and it obviously limits the ability to constrain the MYC model. Nevertheless, fluxes within the MYC model, which describe mycorrhizal C and N exchange with host plants and soil, are reasonably robust against changes in parameter values under unperturbed conditions (Figs. 4, 5, 6). This provides some justification to circumvent the explicit representation of mycorrhizae and instead rely on the use of plant root measurements, i.e. plant fine-root biomass or plant C allocation into fine-root growth, as proxy for mycorrhizal biomass or C transfer from the host plant to fungi. However, the resulting systematic underestimation of simulated flux responses to eCO₂ clearly points to the need to represent mycorrhizae and, thus, also observational data on mycorrhizal fungi. Our study reveals that the organic N uptake capacity ($v_{max}^{\mathrm{m}}$) needs to be better constrained, in particular for saprotrophic mycorrhizae. Even small variations in $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$ affect the ability of mycorrhizae to export N to host plants, which feeds back to plant C exudation (Fig. 3, Eq. A3). Consequently, the MYC model is highly sensitive to $v_{max,N_{\mathrm{SOM}}}^{\mathrm{m}}$ when modelling saprotrophic mycorrhizae. The most important unknown parameters in the representation of decomposing mycorrhizae are the mycorrhiza-induced decomposition acceleration and the general decelerated soil turnover ($f_{{\mathit{\tau }}_{\mathrm{SOM}}}$ and ${\mathit{\tau }}_{\mathrm{SOM}}^{*}$). As these factors do not affect plant–mycorrhiza interactions directly, the MYC model itself is less sensitive to variations, which leads to small variations among LHS simulations, but the combination strongly affects the ecosystem dynamics by potentially degrading almost all SOM.

4.2.2 MYC model limitations due to model constraints

In its current version, the MYC model does not consider P and water supply by mycorrhizal fungi, and QUINCY neglects mycorrhizal biomass in plant leaf-to-root ratio calculations, which not only ensures a certain ratio of C to nutrient acquisition but also plant physical stability. Consequently, QUINCY plants grow roots in order to fulfil physiological and physical requirements that also take up N and make mycorrhizae potentially redundant for plant N acquisition. Therefore, the inclusion of mycorrhizal fungi to P and water supply and to plant leaf-to-root ratio calculations would potentially intensify the meaning of mycorrhizal fungi for host plants and probably reduce uncertainty, especially under changing environmental conditions.

Furthermore, the MYC model assumes that mycorrhizal fungi are independent of the PFT and growth state of the host plants, as long as the plants exude sufficient amounts of C, which does not explicitly vary across PFTs in the current form of the model. This is, due to our approach to include mycorrhizal functionalities and test those, instead of including mycorrhizal types, which one could link to PFTs. Having said that, the plant type and growth state can indirectly affect mycorrhizal function via the extent of N limitation and C available for exudation.

The representation of the SOM as a pool with a homogeneous C : N ratio and chemical properties in QUINCY is a conceptual problem in most state-of-the-art TBMs. In reality, SOM is not a homogeneous pool; rather, it is a complex mixture of substrates with varying C : N ratios, microbial accessibility and decomposability. Describing this variability with only two pools is a major limitation and is particularly problematic for modelling SOM-mining mycorrhizae. As the C : N ratio of the slow SOM pool has to be maintained in the model, uptake of organic N requires a loss of SOC, either through respiration or via incorporation into mycorrhizal biomass through the carbon contained in amino acids as well as other substances. Even if saprotrophic behaviour is occasionally observed (Wu et al., 2005; Treseder et al., 2007), most mycorrhizae solely rely on plant C supply and lack saprotrophic abilities (Lindahl and Tunlid, 2015; Frey, 2019). Employing a more sophisticated representation of soil organic material and its composition and dynamics using more physically based pool structures, stabilisation mechanisms and the dependence of decomposition on microbial activity would be a pathway towards a more accurate and mechanistic modelling of organic N uptake from soils (Wutzler et al., 2017; Kyker-Snowman et al., 2020; Wang et al., 2020; Yu et al., 2020).

4.3 Outlook: possible incorporation of mycorrhizal fungi into TBMs

An increasing number of studies show that mycorrhizal fungi are an important factor for understanding vegetation–soil interactions and soil dynamics as well as for predicting plant responses to elevated CO₂ (e.g. Cheng et al., 2012; Terrer et al., 2018; Shi et al., 2019; Jiang et al., 2020; Braghiere et al., 2022). Our study underlines the importance of considering mycorrhizal fungi in ecosystem models to explain the observed responses of plant N acquisition under eCO₂. Additionally, we show the importance of active plant C exudation to the rhizosphere in order to liberate N from SOM, which allows us to successfully reproduce observed changes in SOC δ¹³C.

Generally, the effects of mycorrhizal fungi on ambient simulations are minor, but the improved eCO₂ simulations suggest that our MYC model can better represent plant control on N acquisition than previous approaches under future conditions. A plant-controlled short-cut of the terrestrial N cycle by accessing SOM as nutrient source is important for reproducing observed responses to eCO₂; thus, this indicates the necessity to incorporate the effects of mycorrhizal fungi on SOM decomposition and plant nutrition into TBMs. However, the addition of any new process into the model always comes with the problems of added complexity, an increased number of parameters and the possibility of equifinality. Particularly given the few available data on mycorrhizae, their inclusion into TBMs can lead to increased uncertainty.

One possible solution would be an implicit formulation of mycorrhizal function in the form of active plant C investment into rhizosphere processes that liberate N from SOM, commonly known as plant priming through root exudates. As we show that the interaction with SOM has the biggest effect on model predictions of the ecosystem response under elevated CO₂, this may in itself improve the representation of plant N acquisition; thus, it has the potential to reduce current uncertainty in modelled future land C uptake. A similar approach has been tested by Fisher et al. (2010), who included organic N uptake in their model.

Another approach to circumvent the problem of too few data could be exploratory modelling studies, where the MYC model with its default parameters is applied globally to qualitatively evaluate the expected effects in scenarios with rising CO₂ (Braghiere et al., 2021). Given the robustness of the MYC model against parameter variations, which we have tested using LHS, we assume that the MYC model would perform well globally. However, we realise that a number of parameters will likely differ with variables such as PFT, soil type and climate; therefore, the quantitative result of such a model experiment would suffer from additional uncertainty given the lack of a suitable source of information for parameters and their global scaling (Braghiere et al., 2022), but the trends and patterns could be a road map for future modelling and data collection.

At the site scale, the inclusion of P acquisition into the MYC model could help to better understand the underlying dynamics and reduce uncertainties, as this would allow us to test the MYC model at P-limited sites or at sites that are both N- and P-limited and, thus, increase the amount of suitable testing sites. We are planning this for future studies.

5 Conclusion

Here, we implemented a novel representation of mycorrhizal fungi into the QUINCY terrestrial biosphere model by representing mycorrhizal functionalities rather that taxonomic differences. We show that, while it does little to improve predictions under ambient conditions, the MYC model improves our predictions under elevated CO₂, most strikingly in the model variant that represents mycorrhizal interactions with soil organic matter. Our study highlights the importance of improving the representation of plant nutrient acquisition strategies in models to improve future predictions of plant growth responses under increased atmospheric CO₂ levels as well as the importance of crucial data needed to parameterise and validate such models.

Appendix A: Detailed model description

A1 Plant C exudation to mycorrhizal fungi

C exudation ($E_{\mathrm{C}}^{\mathrm{p}\mathrm{2}\mathrm{m}}$) from the host plant to mycorrhizal fungi is constrained by the following, as per the concepts of MYCOFON (Meyer et al., 2010):

• the amount of available C, which is freshly assimilated C through photosynthesis (Eq. A2);

• plant N demand, where increasing N demand increases C allocation to mycorrhizal growth (Eq. A3a);

• an interaction term that allows for a decline in exudation if plants do not receive nutrients in return, which could happen either in the case of satisfied demand by plants or in the case of severe N limitation, thereby turning mycorrhizal fungi into competitors that do not deliver any N to host plants;

• a minimum and maximum amount of mycorrhizae compared with plant fine roots, to avoid complete mycorrhizal death and to prevent plants overspending C.

C exudation can be defined as follows:

$$\begin{array}{}\text{(A1)}& E_{\mathrm{C}}^{\mathrm{p}\mathrm{2}\mathrm{m}}=\mathrm{MIN}\left(f_{\mathrm{m}\mathrm{2}\mathrm{r},\min }\times C_{\mathrm{fr}}^{\mathrm{p}}-C^{\mathrm{m}},E_{{\mathrm{C}}_{max}}^{\mathrm{p}\mathrm{2}\mathrm{m}}\right),\end{array}$$

where f_m2r,min is a parameter for the minimal amount of mycorrhizae compared with plant fine roots, $C_{\mathrm{fr}}^{\mathrm{p}}$ is plant fine-root biomass, C^m is mycorrhizal biomass and $E_{{\mathrm{C}}_{max}}^{\mathrm{p}\mathrm{2}\mathrm{m}}$ is the maximum exudation.

We determine the maximum exudation rate using available C as well as the optimal ratio between mycorrhizae and fine roots under current conditions.

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{E}_{{\mathrm{C}}_{max}}^{\mathrm{p}\mathrm{2}\mathrm{m}}& =\mathrm{MIN}(f_{E_{{\mathrm{C}}_{max}}^{\mathrm{p}\mathrm{2}\mathrm{m}}}\times C_{\mathrm{labile}}^{\mathrm{p}}\times {\mathit{\tau }}_{\mathrm{labile}},f_{\mathrm{m}\mathrm{2}\mathrm{r},\mathrm{opt}}\\ & \times C_{\mathrm{fr}}^{\mathrm{p}}-C^{\mathrm{m}}),\end{array}\end{array}$$

where $f_{E_{{\mathrm{C}}_{max}^{\mathrm{p}\mathrm{2}\mathrm{m}}}}$ is the maximum share of available C that plants can allocate to mycorrhizal fungi, $C_{\mathrm{labile}}^{\mathrm{p}}\times {\mathit{\tau }}_{\mathrm{labile}}$ describes the available C, and f_m2r,opt is the optimal ratio between mycorrhizal biomass and plant fine-root biomass.

We restrict the optimal ratio between mycorrhizal biomass and fine-root biomass using a minimal and a maximal ratio and then scale it by plant N demand and mycorrhizal N support. Plant N demand and mycorrhizal N support are used to represent current environmental conditions, which drive plant–mycorrhiza interactions, as plant N demand depends on the plant C uptake by photosynthesis, which is linked to stand-related conditions, such as plant type, successional stage, growing season and atmospheric conditions, and on N availability, which is linked to soil properties and is also depicted by mycorrhizal N support.

$$\begin{array}{}\text{(A3a)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{m}\mathrm{2}\mathrm{r},\mathrm{opt}}& =f_{\mathrm{m}\mathrm{2}\mathrm{r},\min }+(f_{\mathrm{m}\mathrm{2}\mathrm{r},\max }-f_{\mathrm{m}\mathrm{2}\mathrm{r},\min })\times {\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}}\\ & \times {\mathit{\eta }}_{\mathrm{N}},\text{with}\end{array}\text{(A3b)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{N}}={\displaystyle \frac{E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}}{E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}+\mathit{ϵ}}}\end{array}$$

Here, f_m2r,min and f_m2r,max are parameters describing the minimal and maximal ratio between mycorrhizal biomass and plant fine-root biomass. ${\mathit{\zeta }}_{\mathrm{N}}^{\mathrm{p}}$ and η_N are scaling factors between zero and one to describe plant N demand and N support by mycorrhizae respectively. Support is derived from mycorrhizal N export ($E_{\mathrm{N}}^{\mathrm{m}\mathrm{2}\mathrm{p}}$) and an infinitesimal small number (ϵ) to stop C exudation, in the case that plants do not benefit from mycorrhizae over a time period of a month, to technically ensure numerical stability.