Implementing belowground controls on nutrient uptake in ELMv2-SPRUCE improves representation of a boreal peatland ecosystem

Boreal peatlands store 13 %–32 % of the global soil carbon (C) stock, a service dependent on plant-mycorrhizal fungi associations. In these nutrient poor systems, ectomycorrhizal and ericoid mycorrhizal fungi supply up to >80 % of the nutrient requirements of their plant hosts, partly with mined nitrogen (N) and phosphorus (P) from soil organic matter that are otherwise inaccessible to plants. Despite the ecological significance, mycorrhizal associations are only represented in a few land surface or ecosystem models. We modify the peatland branch of version 2 of the Energy Exascale Earth System Land Model (ELMv2-SPRUCE) to replace the default photosynthesis-driven inorganic N and P (NP) uptake process with a more realistic representation of the process via three pathways: (1) direct inorganic NP uptake by uncolonized fine roots, (2) indirect inorganic NP acquisition and (3) indirect NP acquisition from organic sources by mycorrhizal roots. We systematically evaluated the performance of the default and modified models with field observations from a whole ecosystem warming and carbon dioxide fertilization experimental site: Spruce and Peatland Responses Under Changing Environment (SPRUCE), in northern Minnesota, USA. The modified model reduces the underestimation of the growth response of shrubs in the default model to warming from 40 %–80 % to 17 %–35 % and reduces the overall relative absolute error on C fluxes from 1.61 to 1.54 in calibration. Improvements on modeled shrub growths and shrub-moss community net ecosystem exchanges are also seen in validation. The improved growth response of shrubs to warming is accompanied by several-fold increase in direct inorganic NP uptake and decrease in fungal colonization rate. The modified model simulates a smaller magnitude of transition of the ecosystem from C sink to C source under warming due to alleviation of plant nutrient limitation. Equifinality analysis shows the newly added parameters in the modified model can be constrained by the observed C fluxes. Sensitivity analysis shows the newly added parameters have stronger statistical interactions than the preexisting parameters in the default model. Overall, the modified model is an improvement over the default ELMv2-SPRUCE and will be a useful tool for understanding boreal peatland change.

1 Background

Boreal peatlands store an estimated 234–546 Gt carbon (C), equal to 13 %–32 % of the global soil C stock (Friedlingstein et al., 2022; Loisel et al., 2017). The high C storage arises from slow decomposition rates driven by the cold, waterlogged, nutrient-limited, and acidic conditions of these ecosystems (Dise, 2009; Frolking et al., 2011; Salmon et al., 2021). Ongoing rapid warming in the northern high latitudes is expected to shift ecosystem C balance, but the magnitude of change remains highly uncertain due to poorly constrained temperature sensitivities of vegetation productivity and soil C decomposition (Ito et al., 2020). More accurate modeling of the mechanisms governing C cycling in boreal peatlands will improve our ability to project future changes in this ecosystem and its feedback to the Earth system.

Among the various biotic and abiotic mechanisms underlying boreal peatland C cycling, plant-mycorrhizal associations represent a key component owing to their central role in nutrient cycling (Shao et al., 2022, 2023b; Shi et al., 2015, 2021). Mycorrhizal fungi deliver nutrients to plants in exchange for carbon and have three broad classes: ectomycorrhizae (EcM), ericoid mycorrhizae (ErM), and arbuscular mycorrhizae (AM). Unlike AM, which are more common in low latitudes and only acquire inorganic nutrients, EcM and ErM can acquire nutrients from soil organic matter (SOM), making them suited to cold, nitrogen (N)-limited ecosystems with slow decomposition rates (Egerton-Warburton et al., 2013; Ward et al., 2022). One estimate suggests that EcM are associated with >75 % of the aboveground plant biomass in the non-permafrost boreal region and ∼ 50 % in the permafrost region; ErM, which selectively colonize ericaceous shrubs, are associated with ∼ 20 % aboveground plant biomass in the permafrost region (Soudzilovskaia et al., 2019). Plants associated with EcM transfer on average ∼ 13 % (0 %–50 % in range) of their net primary productivity (NPP) to fungal symbionts while plants associated with ErM transfer on average ∼ 3.5 % (0 %–14 % in range; Hawkins et al., 2023). The fraction of plant N supplied by EcM or ErM in return varies from <30 % to >80 %, depending on site and plant species, though less is known about the fraction of EcM- or ErM-supplied phosphorus (P; Hilman et al., 2024; Hobbie and Hobbie, 2006; Yin et al., 2022). Beyond nutrient supply, mycorrhizal fungi regulate SOM turnover by competing with free-living saprotrophs, transporting C away from the rhizosphere, and promoting soil aggregate formation and stabilization (Fernandez and Kennedy, 2016; Hawkins et al., 2023; Smith and Wan, 2019).

To date, only a limited number of land surface models – here defined as the land component of Earth system models – simulate mycorrhizal associations (Warren et al., 2015). They generally focus on AM and EcM associations and use the return-on-investment principle, where “return” refers to gains in growth or nutrient uptake, and “investment” refers to the C costs of acquiring nutrient through different pathways (Brzostek et al., 2014). For example, the Community Land Model (CLM) and Energy Exascale Earth System Land Model (ELM) have been linked with the Fixation and Uptake of Nitrogen (FUN) model (Braghiere et al., 2022; Brzostek et al., 2014; Shi et al., 2016). Plants minimize their C expenditure on N and P (NP) uptake by optimally allocating C to satisfy their NP demands among biological fixation, retranslocation, nonmycorrhizal passive and active uptake, EcM uptake, and AM uptake, each of which has a unique C cost function (Braghiere et al., 2022; Brzostek et al., 2014; Shi et al., 2016). Simulations of ELM-FUN suggest that the EcM and AM pathways together supply ∼ 75 % of plant N and ∼ 41 % of plant P globally, and account for ∼ 50 % of the NP uptake-related C costs, but neither ELM-FUN or CLM-FUN consider organic nutrient mining (Braghiere et al., 2022; Shi et al., 2016). The Symbiotic Nitrogen Acquisition by Plants (SNAP) model, which is linked to the Geophysical Fluid Dynamics Laboratory land model LM3 (GFDL-LM3), improves FUN by dynamically simulating fungal biomass, fungal organic nutrient mining, and the resulting C cost to the plants (Sulman et al., 2019). Simulations of GFDL-LM3-SNAP show that EcM-mining of organic N explains the stronger positive response to carbon dioxide (CO₂) fertilization in EcM-dominated ecosystems than AM-dominated ecosystems (Sulman et al., 2019). Also, allowing plants to shift in N uptake pathways results in four times the terrestrial C sequestration relative to fixed N uptake pathways under a 100 ppm increase in atmospheric CO₂ concentration (Sulman et al., 2019).

Terrestrial ecosystem models not coupled to Earth system models have represented mycorrhizal associations in more detail than the return-on-investment models described above. For example, the McGill Wetland Model (MWM) focuses on interactions among moss, ericaceous shrub, and ErM in peatland ecosystems and shows that the shrub-ErM association explain the increased shrub growth and decreased moss growth in a NP fertilization experiment (Shao et al., 2022, 2023b). The MWM explicitly models microbial and ErM biomass dynamics and ErM mining of organic nutrients. The MWM models the shrub-ErM interactions as “excess fluxes” in which (1) shrub transfers C to ErM when shrub C reserve exceeds a set fraction of its total stem and root C, (2) ErM fungi transfers NP to the shrub when the NP contents of the ErM exceed predefined fractions of the C content of the ErM (Shao et al., 2023b). Compared to return-on-investment, the excess flux mechanism may better describe EcM and ErM exchanges with the host plant, because the reciprocity of EcM and ErM are more strongly affected by environmental, developmental, and physiological factors than AM (Bunn et al., 2024; Garcia et al., 2015). The drawback is a large number of parameters. The CoupModel has been used to compare three representations of nutrient limitation: fixed limitation, implicit EcM, and explicit EcM across a climate and N-deposition gradient of EcM-dominated boreal forests (He et al., 2018). In the fixed limitation approach, plant growth is scaled down by a constant nutrient limitation factor throughout the year. The implicit approach omits the EcM intermediary, simulating plant acquisition of NP from organic sources as a function of soil organic nutrient content, plant demand, and optionally root distributions (He et al., 2018, 2021; Svensson et al., 2008). The explicit approach simulates EcM biomass dynamics and organic nutrient mining, with plant transfer of C to EcM determined by belowground allocation, and EcM transfer of NP to plant co-determined by plant demand and excess flux (He et al., 2018, 2021). The implicit and explicit approaches outperform the fixed N limitation approach, and both indicate declining plant dependence on organic N from the more N-limited northern Sweden to the less N-limited southern Sweden; yet the explicit parameterization is more difficult to constrain, and the implicit and explicit approaches differ in the simulated litter production, soil respiration, and the magnitude of the north-south trend (He et al., 2018).

The above reviewed modeling studies demonstrate that mycorrhizal associations are needed for more accurate simulation of nutrient limitation on productivity and the resulting feedback to Earth system and land surface models can benefit from testing alternative model structures than return-on-investment schemes and understanding the parameterization difficulty. We address this research gap by adding implicit representation of EcM and ErM associations into the NP uptake processes of a peatland branch of ELM, ELMv2-SPRUCE (Griffiths et al., 2017; Shi et al., 2015, 2021). AM is not added because it is not a key component of northern peatland ecosystems (Egerton-Warburton et al., 2013). We compare the original model, hereafter “ELM-OLD”, and the modified model, hereafter “ELM-MYCI” (for mycorrhizal-implicit), against observed C fluxes, pore water NP concentrations, resin-exchange-measured plant available NP, and peat C-N-P stocks from the Spruce and Peatland Responses Under Changing Environment experiment (SPRUCE; Griffiths et al., 2017; Griffiths and Sebestyen, 2016; Hanson et al., 2020a; Iversen et al., 2022; Salmon et al., 2021). SPRUCE is a whole ecosystem warming and CO₂ fertilization experiment located in a boreal peatland ecosystem (Hanson et al., 2017). The site has EcM-associated trees (black spruce [Picea mariana] and tamarack [Larix laricina]) and various species of ErM-associated ericaceous shrubs, offering an array of interactions between plants, fungi, and soils under experimental treatments that have not been tested by the above-reviewed modeling studies. Warming at SPRUCE increased shrub productivity, decreased Sphagnum moss productivity, and increased resin-exchange nutrient availability (Hanson et al., 2020a, 2025; Iversen et al., 2022). A persistent issue in ELM-OLD has been the inability to reproduce the larger increase in shrub productivity relative to tree productivity under warming (Shi et al., 2021). Using ELM-MYCI, we test the hypothesis that the shrub responses can be explained by decreasing dependence on ErM in response to higher nutrients availability under warming, akin to the findings or suggestions of multiple previous studies (Defrenne et al., 2021; Duchesneau et al., 2024a; He et al., 2018; Shao et al., 2023b).

2 Data and Methods

2.1 Site Description

The SPRUCE experimental site is within the S1 Bog of the United States Department of Agriculture Forest Service Marcell Experimental Forest, located in northern Minnesota, USA (47°30.476^′ N, 93°27.162^′ W, 418 m above mean sea level; Hanson et al., 2017, 2020a; Kolka et al., 2011; Salmon et al., 2021). The bog is an acidic, raised-dome ombrotrophic bog with nutrient inputs only from atmospheric deposition and nitrogen (N) fixation. The open forest canopy at the site is Picea mariana (Mill.) B.S.P. (black spruce) with occasional Larix laricina (Du Roi) K. Koch (eastern tamarack). The trees were harvested in strip cuts in 1969 and 1974 and the current canopy is mostly regenerated from the 1974 strip cut (Hanson et al., 2016a). The understory is dominated by ericaceous shrubs (Rhododendron groenlandicum [Oeder] Kron & Judd [Labrador tea], Chamaedaphne calyculata [L.] Moench. [leatherleaf], Vaccinium angustifolium Aiton [lowbush blueberry], Vaccinium oxycoccos L [bog cranberry], and Kalmia polifolia [bog laurel]) with a small biomass pool of forbs and sedges (Hanson et al., 2025). The bryophyte layer is dominated by Sphagnum spp. This vegetation community is represented by the following plant functional types (PFTs) in both ELM-OLD and ELM-MYCI: boreal evergreen needleleaf for black spruce, boreal deciduous needleleaf for tamarack, boreal deciduous shrub for the ericaceous shrubs, and Sphagnum moss. The forbs and sedges are not modeled but comprise less than 10 % understory cover according to pretreatment surveys (Iversen et al., 2017b).

The detailed whole ecosystem warming and CO₂ fertilization experimental setup is reported elsewhere (Hanson et al., 2017). Briefly, the experiment has two unenclosed, ambient plots (ambient temperature, ambient CO₂) and five pairs of enclosures that target five whole ecosystem warming levels (+0, +2.25, +4.5, +6.75, and +9 °C) above ambient temperatures crossed with ambient and elevated (+500 ppm) CO₂. The belowground heating extends 3 m into the peat profile and began in June 2014. The aboveground warming began in August 2015. The CO₂ fumigation began on 15 June 2016.

2.2 The default ELMv2-SPRUCE model (ELM-OLD)

ELM is the land component of the Energy Exascale Earth System Model (E3SM), which consists of atmosphere, land, ocean, sea ice, and land ice components (Burrows et al., 2020; Yang et al., 2019, 2023). ELM-OLD is currently branched off ELM version 2 with improved peatland processes, including hummock-hollow hydrological interactions (Shi et al., 2015) and the Sphagnum moss PFT (Shi et al., 2021). ELM-OLD has been used primarily for site-level simulations, in which we represent the bog as two interacting grid cells that represent a hummock soil column and a hollow soil column (Shi et al., 2015). Each soil column has multiple PFTs that compete for water and nutrients (Shi et al., 2021). Soil decomposition uses a first-order decay model with one coarse woody debris pool, three plant litter pools (labile, cellulose, lignin), and four SOM pools (Burrows et al., 2020; Oleson et al., 2013; Supplement, Sect. S1.1.3). Belowground nutrient competition among the PFTs and the soil decomposition process is simulated with the Relative Demand approach (Sect. S1.1.2; Fig. 1a). Plant photosynthesis creates potential NP uptake through the fixed C : N and C : P in the plant structural tissues (leaf, fine root, coarse root, stem). Soil decomposition creates potential NP immobilization due to the need for extra NP when C decomposes from an upstream pool that has higher C : N and C : P into a downstream pool that has lower C : N to C : P. The potential uptakes and immobilization are compared to the total available inorganic NP in soil, and scaled down by the same factor so that the total available inorganic NP is not exceeded (Burrows et al., 2020; Thornton and Rosenbloom, 2005). Due to this NP limitation, some photosynthesized C cannot become growth in the structural tissues; the extra C enters the nonstructural carbohydrates (NSC) pool as C reserve (Burrows et al., 2020).

Figure 1 Simplified illustration of the belowground nutrient competition between the vascular plant functional types and the soil decomposition processes in ELM-OLD (a) and ELM-MYCI (b). Boxed arrows indicate potential plant uptakes. Dashed arrows indicate directions of nutrients flow, with blue for mineral nutrients and purple for nutrients locked in organic forms. Solid arrows between the decomposition pools indicate direction of carbon flow, with unattached arrow heads meaning respirational losses. CWD – coarse woody debris.

2.3 The modified ELMv2-SPRUCE model (ELM-MYCI)

ELM-MYCI is designed to improve the process-realism of nutrient uptake for the three vascular PFTs (spruce, tamarack, shrubs) by considering the following three pathways: (1) direct inorganic nutrient uptake by uncolonized fine roots, PATH^root, (2) indirect inorganic nutrient acquisition by mycorrhizal roots, PATH^myc,inorg, and (3) indirect nutrient acquisition from organic sources by mycorrhizal roots, PATH^myc,org (Fig. 1b). The entire set of equations and detailed descriptions are provided in Sect. S1.1.4–S1.1.10. The major equations and assumptions are described here. Unless otherwise noted, the equations are only shown for N, and the equations for P can be obtained simply by replacing all N with P.

Like in ELM-OLD, ELM-MYCI still calculates the NP demand implied by photosynthesis, but decouples this demand from the potential NP take, which are determined through the three pathways. ELM-MYCI uses fungi-colonization fraction to idealize the fine root into a uncolonized part, which can only use PATH^root, and a colonized part, which can only use PATH^myc,inorg and PATH^myc,org. The colonization fraction of EcM for the two tree PFTs and ErM for the shrub PFT is a clipped linear function of annual average soil inorganic N (Eq. 1). The sensitivity to P is currently set to zero in Eq. (1), because observational constraints on this sensitivity are weak in past experimental studies (Bashian-Victoroff et al., 2025) and at the SPRUCE site (Sect. S1.1.5).

The potential inorganic NP uptake via PATH^root is dependent on uncolonized fine root surface area, soil inorganic NP concentrations, soil temperature, soil moisture, and the current NP-limitation level of the plant (Eq. 2). The fine root surface area, $A_{\mathrm{froot},i,j}$, is calculated from modelled fine root biomass, and observed fine root radius and density at the SPRUCE site (Eq. S32 in the Supplement). The dependence on soil inorganic NP concentrations, F_j(N_conc,i), follows the conventional Michaelis-Menten form (Eq. S22; Knox et al., 2024). The dependence on soil temperature, F(T_soi,i), follows the conventional Q₁₀ form (Eq. S13). The dependence on soil moisture, F(Θ_soi,i), follows a previous unimodal function to let both dry and excessively wet soil inhibit nutrient uptake (Eq. S14; Frolking et al., 2002). The NP-limitation feedback term, F(F_Nlimit,j), is exactly one when the PFT is not NP-limited and decreases asymptotically to zero as soil inorganic NP become more abundant (Eq. S15; Fig. S3 in the Supplement). This term prevents numerical instability in the model and is somewhat supported by experimental observations (Glass et al., 2002).

The potential inorganic NP acquisition via PATH^myc,inorg is dependent on colonized fine root biomass, the same soil inorganic NP concentrations, soil temperature, and soil moisture multipliers as PATH^root, and NSC availability (Eq. 3). The NSC availability term reflects fungal dependence on C transfer from the host plant (Eq. S16; Bunn et al., 2024; He et al., 2018; Shao et al., 2023b). The potential NP acquisition from organic sources via PATH^myc,org is dependent on colonized fine root biomass, the same soil temperature and moisture multipliers as PATH^root and PATH^myc,inorg, and the same NSC availability term as PATH^myc,inorg (Eq. 4). PATH^myc,org does not include a Michaelis-Menten term for organic NP concentrations, because the peatland environment is rich in organic substrates. Also, the soil decomposition pools in ELM are rate-based rather than corresponding to actual chemical compounds and do not distinguish between solid and dissolved phases (Oleson et al., 2013), making it difficult to separate fungi-accessible from fungi-inaccessible fractions (Näsholm et al., 2009; Talbot and Treseder, 2010). The total potential mycorrhizal inorganic and organic NP acquisitions is capped by 50 % of net photosynthesis per time step, scaled by a constant C-cost factor (Eq. S25). The 50 % threshold was selected based on the maximum percentage of C allocated to mycorrhizal fungi in a previous meta-analysis (Hawkins et al., 2023). The cap, like the NSC availability term, reflects fungal dependence on C input from the host plants.

The potential rates of the three pathways are compared to soil inorganic and organic NP availability to determine actual rates. The potential inorganic NP uptake or acquisition from PATH^root, PATH^myc,inorg, and soil decomposition processes are scaled down by soil inorganic NP availability analogously to the default Relative Demand approach (Eq. S3, S4, S9, S11). The potential organic NP acquisition from PATH^myc,org per hourly time step is capped by 0.0001 of the total organic NP in the labile, cellulose, and lignin plant litter pools (Fig. 1; Eq. S19–S20). The 0.0001 multiplier is adopted from the upper bound in CoupModel (He et al., 2018), and is intentionally chosen to be a large number that prevents unrealistically large uptake rate rather than serving as a true upper bound on the fungi-accessible fraction of organic NP. The restriction to only access plant litter pools is again due to structural constraint in ELM's soil decomposition pools. The plant litter pools allow variable C : N and C : P ratios, but the SOM pools do not, so that mycorrhizal acquisition of N and P from the SOM pools will require considering how much C is lost as fungal respiration, which is beyond the scope of this project (Sect. S1.1.7). Finally, the actual rates are compared to the implied demand by photosynthesis to determine plant structural growth. The actual rate of PATH^myc,org is used to subtract NP from the three plant litter pools, distributing the acquisition proportional to the size of each pool.

$$\begin{array}{}\text{(1)}& {\displaystyle }{M}_{\mathrm{myc},j}=max\left(\mathrm{0},min\left(\mathrm{1},a_j+b_j{N}_{\mathrm{soil},\mathrm{annavg}}\right)\right)\text{(2)}& {\displaystyle }\begin{array}{rl}{N}_{\mathrm{froot},i,j}& =v_{\mathrm{N},\mathrm{froot},j}\left(\mathrm{1}-M_{\mathrm{myc},j}\right)A_{\mathrm{froot},i,j}{F}_j\left(N_{\mathrm{conc},i}\right)\\ & \cdot F\left(T_{\mathrm{soi},i}\right)F\left({\mathrm{\Theta }}_{\mathrm{soi},i}\right)F\left(F_{\mathrm{Nlimit},j}\right)\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{N}_{\mathrm{myc},\mathrm{pot},\mathrm{inorg},i,j}& =v_{\mathrm{N},\mathrm{myc},j}{M}_{\mathrm{myc},j}{C}_{\mathrm{froot},j}{F}_{\mathrm{froot},i,j}\\ & \cdot F_j\left(N_{\mathrm{conc},i}\right)F\left(T_{\mathrm{soi},i}\right)F\left({\mathrm{\Theta }}_{\mathrm{soi},i}\right)\\ & \cdot F\left(F_{\mathrm{Nlimit},j}\right)F\left(C_{\mathrm{ns},j}\right)\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}{N}_{\mathrm{myc},\mathrm{pot},\mathrm{org},i,j}& =u_{\mathrm{N},\mathrm{myc},j}{M}_{\mathrm{myc},j}{C}_{\mathrm{froot},j}{F}_{\mathrm{froot},i,j}\\ & \cdot F\left(T_{\mathrm{soi},i}\right)F\left({\mathrm{\Theta }}_{\mathrm{soi},i}\right)F\left(F_{\mathrm{Nlimit},j}\right)\\ & \cdot F\left(C_{\mathrm{ns},j}\right)\end{array}\end{array}$$

•  

  i – soil layer index

•  

  j – PFT index

•  

  M_myc,j – fraction of fine roots colonized by EcM (for the spruce and tamarack PFTs) or ErM (for the shrub PFT)

•  

  N_soil,annavg – annual average soil inorganic N (${\mathrm{NH}}_{\mathrm{4}}^{+}+{\mathrm{NO}}_{\mathrm{3}}^{-}$) content in the rooting zone, g N m⁻³

•  

  a_j – intercept parameter

•  

  b_j – slope parameter

•  

  $N_{\mathrm{froot},i,j}$ – the potential inorganic N uptake rate via PATH^root, g N m⁻² ground area s⁻¹

•  

  $N_{\mathrm{myc},\mathrm{pot},\mathrm{inorg},i,j}$ – the potential inorganic N acquisition rate via PATH^myc,inorg, g N m⁻² ground area s⁻¹

•  

  $N_{\mathrm{myc},\mathrm{pot},\mathrm{org},i,j}$ – the potential N acquisition rate from organic sources via PATH^myc,org, g N m⁻² ground area s⁻¹

•  

  $v_{\mathrm{N},\mathrm{froot},j}$ – the maximum inorganic N acquisition rate per unit uncolonized fine-root surface area, g N m⁻² ground area s⁻¹

•  

  $v_{\mathrm{N},\mathrm{myc},j}$ – the maximum inorganic N acquisition rate per unit colonized fine-root biomass, g N g C⁻¹ s⁻¹

•  

  $u_{\mathrm{N},\mathrm{myc},j}$ – the maximum organic N acquisition rate per unit colonized fine-root biomass, gN g C⁻¹ s⁻¹

•  

  $A_{\mathrm{froot},i,j}$ – total fine root surface area in one soil layer, cm² m⁻² ground area

•  

  C_froot,j – total fine root biomass in one soil layer, g C m⁻² ground area

•  

  F_j(N_conc,i) – Michaelis-Menten multiplier of soil inorganic N concentration (N_conc,i, g N m⁻³ soil volume)

•  

  F(T_soi,i) – Q₁₀ multiplier of soil temperature (T_soi,i, °C)

•  

  F(F_Nlimit,j) – a feedback factor to prevent infinite N uptake when inorganic N is abundant (F_Nlimit,j is the PFT's N-limitation level in the previous time step)

•  

  F(C_ns,j) – degree of NSC saturation in the plant (C_ns,j is the NSC biomass in the plant, g C m⁻² ground area)

2.4 Simulation protocol

Following previously established protocols (Griffiths et al., 2017; Hanson et al., 2020a; Shi et al., 2021), we first conducted a single-grid simulation that consists of an accelerated spin-up of 207 years, a normal spin-up of 407 years, and an 1850–2014 transient simulation, and then branched the simulations into eleven treatments corresponding to one control simulation for unenclosed plot + five pairs of enclosures (Sect. 2.1) during 2015–2023. The control simulation only uses one of the two unenclosed plots, which is labeled plot 7 in the experiment (Hanson et al., 2020a), because it has a longer water table record, which is needed to force ELMv2-SPRUCE. The accelerated and normal spin-ups were driven by cyclic ambient meteorological forcing during 2015–2023, preindustrial CO₂ concentration, preindustrial N deposition, and constant land cover. The transient simulation cyclically used the ambient meteorological forcing during 2015–2023, historically varying CO₂ concentration and N deposition, and included the 1974 strip cut event where 99 % aboveground tree biomass was removed. The treatment simulations were forced by meteorological observations and water table depths in each plot during 2015–2023 (Hanson et al., 2016b, 2020b). The simulated water table depths in the two columns equilibrate with each other and observed water table depths (Shi et al., 2015). The atmospheric CO₂ concentrations in the elevated CO₂ plots were set to 500 ppm above ambient level starting from 15 March 2016. Within the grid, the hummock soil column was set to 64 % of the area and hollow 36 % (Graham et al., 2020). A limitation of this version of ELM is that we do not represent multiple canopy layers. Therefore, we must specify fractional coverages for each PFT that add to 100 % total. We started with the default assumption that each of the four PFT covers 25 % and then adjusted for the observed distribution of the two tree types. Within each soil column, the PFT fractions were: needleleaf evergreen boreal tree 36 % (for spruce), needleleaf deciduous boreal tree 14 % (for larch), broadleaf boreal deciduous shrub 25 % (for ericaceous shrubs), and Sphagnum moss 25 % in the pre-treatment simulations and fractionally adjusted using the annually observed fractional coverages in the treatment simulations (Table S8 in the Supplement). All the simulations used an hourly time step.

2.5 Model calibration data

For the parameter optimization experiments described in Sect. 2.7.1, we focused on the following annual C fluxes: (1) the aboveground NPP of the spruce trees (AGNPP_spruce), (2) the aboveground NPP of the tamarack trees (AGNPP_tamarack), (3) the aboveground NPP of the shrubs (AGNPP_shrub), (4) the NPP of Sphagnum moss (NPP_moss), (5) the belowground NPP of total tree and shrub fine roots (BGNPP_treeshrub), (6) heterotrophic respiration (HR; Hanson et al., 2018a, b, 2020a; Norby et al., 2019). We additionally summed up (1–5) to obtain (7) an aggregated “NPP” term – note this is not true NPP because it does not contain coarse root production. There are some temporal mismatches across the datasets: BGNPP_treeshrub observations only span 2016–2017 (Malhotra et al., 2020b); the other observations span 2016–2021, but year 2020 was excluded due to the high uncertainty associated with the limited measurements taken during the COVID era (Hanson et al., 2020a; Norby and Childs, 2018). To facilitate concise comparison, we summarized each of those variables into two mean values and two temperature sensitivities, similar to a previous approach at the site (Hanson et al., 2020a). The mean values were calculated, respectively, over all the years in the ambient CO₂ plots and over all the years in the elevated CO₂ plots. The temperature sensitivities were calculated as the slope of least squares linear regression between each C-flux variable and observed mean annual 2 m air temperatures, respectively, over all the years in the ambient CO₂ plots, and over all the years in the elevated CO₂ plots. Those means and slopes were used in parameter optimization objective function (Eq. 5).

2.6 Model evaluation data

The datasets described in the following subsections were withheld from model calibration and used solely for evaluation.

2.6.1 Post-calibration annual C fluxes

We reserved the following annual C fluxes for post-calibration model evaluation: AGNPP_spruce and AGNPP_tamarack, available for year 2022 in the ten treatment enclosures; AGNPP_shrub, available for year 2023 in all eleven plots; the growing season average NEE of the shrub-moss community (NEE_shrubmoss, gC m⁻² gs⁻¹, “gs”), available for year 2023 in the ten treatment enclosures. The data sources of AGNPP_spruce, AGNPP_tamarack, and AGNPP_shrub are the same as in Sect. 2.5. The growing season mean NEE_shrubmoss was calculated from the average of gapfilled 15 min automated chamber measurements made on the shrub-moss community (Stelling et al., 2024). Growing season is defined as 1 May to 31 October, to match the observational period (27 April to 27 October 2023). The gapfilling method is adapted from previous work at the SPRUCE site (Walker et al., 2017) and is described in detail in Sect. S1.3. Note the automated chamber measurements and the gap-filling method used here are preliminary and were applied solely for the purposes of this study, pending availability of the full analysis of the automated chamber data (J. Stelling, personal communication, 2006). To obtain the modelled equivalence of NEE_shrubmoss, we considered what component fluxes were included in the observed NEE_shrubmoss. The chambers only enclosed shrubs and moss aboveground. Therefore, the chamber fluxes included aboveground shrub and moss GPP and respiration. The belowground collars of the chambers did not exclude tree roots. Therefore, the chamber fluxes included belowground soil HR, tree root respiration, shrub root respiration, and belowground moss respiration. Furthermore, because the chambers have small footprints, the sampled shrub and moss growths were often not representative of the average plot conditions. Therefore, we down-scaled the modelled plot-level shrub and moss GPP and respiration to the chamber-level using each PFT's ratio of measured aboveground biomass in the chamber (gC m⁻²) to the measured aboveground biomass in the corresponding plot (gC m⁻², Table S8). We assumed the soil HR and tree root respiration components of the chamber fluxes correspond to plot-level averages due to lack of prior information. In summary, the modelled equivalence of NEE_shrubmoss is calculated as (Plot-average HR + Plot-average tree root respiration + Scaled shrub autotrophic respiration + Scaled moss autotrophic respiration – Scaled shrub GPP – Scaled moss GPP).

2.6.2 Annual maximum leaf area index

To evaluate the simulated C biomass, we compared the models against the annual maximum leaf area index (LAI) of the two trees and shrub, measured using LICOR LAI 2200 device during 2015–2020 (McPartland et al., 2019). The annual C fluxes and LAI observations all have direct correspondence with modeled variables. AGNPP_spruce and AGNPP_tamarack have strong pre-treatment variation that impacts the interpretation of results (Hanson et al., 2025). Therefore, we fitted ordinary least squares linear regression models to remove the pre-treatment effect; because LAI is closely related to aboveground NPP, we used the same method to remove the pre-treatment variation in the LAI of the two tree species (Sect. S1.2). We did not compare the models against observed net C exchange or methane flux data because the net C exchange has small components of lateral outflow of organic C and dissolved inorganic C (Hanson et al., 2020a), which are not yet considered in ELMv2-SPRUCE, and the model cannot yet fully capture methane dynamics.

2.6.3 Porewater nutrient concentrations

For nutrient cycling, we evaluated the models against several forms of observations. The first of those was porewater ammonia (NH${}_{\mathrm{4}}^{+}$), nitrate (NO${}_{\mathrm{3}}^{-}$), and inorganic P (PO${}_{\mathrm{4}}^{\mathrm{3}-}$) concentrations, which characterize the pool of dissolved nutrients in soil. The porewater concentrations were measured roughly twice per month during the non-frozen months of 2015–2020 near the bottom of the rooting profile (30 cm) in the hollow (Griffiths et al., 2016). The pore water nutrient concentrations were measured on a per water volume basis (mg N L⁻¹, mg P L⁻¹), but ELM can only simulate nutrient concentrations on a per soil volume basis (gN m⁻³ soil, gP m⁻³ soil). Therefore, we divided the simulated NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, and soluble P (as defined in Yang et al., 2019) concentrations by the simulated volumetric soil water content to bring them to the same unit as the observations. Because the sampling interval was irregular and relatively infrequent, and nutrient concentrations in water can exhibit sharp, episodic fluctuations (Basu et al., 2010), we chose not to interpolate the observations to annual level like the annual C fluxes (Hanson et al., 2018a, b, 2020a; Norby et al., 2019). Instead, we down sampled the model outputs to only include the dates on which pore water concentrations were measured, ensuring that both observed and modeled data reflect the same seasonality. For each enclosure, we then summarized the observed and down-sampled modeled porewater nutrient concentrations into two statistics: the mean value and the temporal linear trend, to facilitate concise comparison. Here, the temporal linear trend is defined as the slope of a least-squares regression between the modeled/observed values and the elapsed days since 1 January 2015. We chose enclosure-specific temporal trends, instead of cross-enclosure temperature sensitivity (as done in Sect. 2.5), because the former metric captures well the behavior of the data when plotted as a time series, whereas scatterplots between the observed porewater nutrient concentrations and the air temperature on the corresponding dates do not show clear correlation.

2.6.4 Resin-exchange plant available nutrients

The second form of nutrient cycling observations we used is resin-exchange NH${}_{\mathrm{4}}^{+}$ and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (Iversen et al., 2017a). Those values approximately represent plant-available nutrients but previous study suggests that in shallow soil layers (10 cm), competition with roots invalidates this interpretation (Iversen et al., 2017a, 2022). Therefore, we focus on the resin-exchange nutrients measured at 30 cm in both the hummock and hollow, and at 60 cm in the hollow during 2015–2018, at which depths few roots exist in the anaerobic peatland environment (Iversen et al., 2017a). The raw measurements were cumulative absorbed weights of NH${}_{\mathrm{4}}^{+}$ or PO${}_{\mathrm{4}}^{\mathrm{3}-}$ per resin capsule surface area during roughly monthly intervals; the values used in this study are the aggregated annual averages by the data collector (Iversen et al., 2017a). Following previous study, we compared the temperature sensitivities of the observed resin-exchange nutrients to the temperature sensitivities of the most comparable modeled variable, annual average net nutrient mineralization rates (NET_NMIN [gN m⁻³ s⁻¹], NET_PMIN [gP m⁻³ s⁻¹]) for a qualitative comparison (Iversen et al., 2022). To remove unit difference between the observed and modeled quantities, we divided the observed values and the modeled values, respectively, by the mean value of each during the observational period, separately for each nutrient species and soil depth. The resulting regression slopes against temperature represent relative (mean-scaled) changes per unit change in soil temperature. This normalization approach is commonly used in elasticity analysis for comparing effects across variables with different units (Sydsaeter and Hammond, 1995).

2.6.5 Pretreatment peat C, N, P stocks

The final observation of nutrient cycling we used for evaluation are pre-treatment vertical profiles of C, N, and P stocks in the peat (Griffiths et al., 2017; Salmon et al., 2021). The compared model variables are the summed C, N, and P contents of all four SOM pools (1, 2, 3, 4) in all the soil layers in the decomposition processes (Sect. S1.1.3; Oleson et al., 2013) in the control plot.

2.7 Sensitivity experiments

2.7.1 Parameter optimization

In both ELM-OLD and ELM-MYCI, the soil decomposition process uses column-level parameters that are shared between hummock and hollow. Vegetation processes like photosynthesis, respiration, and nutrients uptake/acquisition use PFT-specific parameters that are also shared between hummock and hollow. We used observed parameter values as much as possible in ELM-OLD (Table S2) and ELM-MYCI (Table S4). For those unobserved, we used either optimized values obtained from 4000-member ensemble simulations that are described below (Tables S3 and 1) or manually set default values (Tables S2 and S4) if a parameter has relatively minor influence on the quantities of interest in this study. All the parameter optimization simulations excluded moss parameters for the purpose of consistency – because moss is not modified in this study, the newly added parameters do not cover moss.

Table 1 Newly added parameters in ELM-MYCI that are optimized.

We selected eleven preexisting parameters in ELM-OLD to optimize based on previous sensitivity findings (Meng et al., 2021; Ricciuto et al., 2018). Those parameters affect photosynthesis, plant respiration, plant allocation, and soil decomposition processes (Table S3). For each parameter, we generated an equal number of uniform random samples between predetermined upper and lower bounds from previous studies (Griffiths et al., 2017; Meng et al., 2021; Ricciuto et al., 2018; Table S3). We then formed all possible combinations of these samples across the parameters to create a total of 4000 samples. We then ran the model with these 4000 parameter sets. Finally, we ranked these parameter sets by relative absolute error (RAE) and selected the sample with the lowest RAE as the optimized parameter values. The RAE formula is shown in Eq. (5) and considers the relative errors in mean and in temperature sensitivity of annual C fluxes:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{RAE}& ={\displaystyle \frac{\mathrm{1}}{\left|V\right|}}\sum _{\begin{array}{l}v\in V\\ c\in \left\{\mathrm{ambient},\mathrm{elevated}\right\}\end{array}}\\ & \cdot \left(\left|{\displaystyle \frac{A_{\mathrm{sim},\mathrm{v},\mathrm{c}}-A_{\mathrm{obs},\mathrm{v},\mathrm{c}}}{U_{\mathrm{obs},\mathrm{v},\mathrm{c}}}}\right|+\left|{\displaystyle \frac{S_{\mathrm{sim},\mathrm{v},\mathrm{c}}-S_{\mathrm{obs},\mathrm{v},\mathrm{c}}}{Q_{\mathrm{obs},\mathrm{v},\mathrm{c}}}}\right|\right)\end{array}\end{array}$$

where the subscripts sim means simulated, obs means observed, v is variable, c is the CO₂ treatment (ambient or elevated), V is the set of variables compared, |V| is the number of variables, A is the mean magnitude of the variable across all the treatments and years, S is the linear regression slope of the variable against annual mean air temperature across all the treatments and years, U is the observational uncertainty in mean, Q is the observational uncertainty in slope. During calibration, the performance metric calculated on the annual C fluxes that are items 1–6 in Sect. 2.5, because those are directly comparable to model outputs. Other observations are either reserved for validation (Sect. 2.6.1–2.6.2), not directly comparable to model outputs (Sect. 2.6.3–2.6.4), or pretreatment (Sect. 2.6.5). We did not include the aggregate NPP (item 7) or annual maximum LAI because of strong overlap with items 1–6. The $U_{\mathrm{obs},\mathrm{v},\mathrm{c}}$ values were estimated from pre-treatment standard deviation across the enclosure locations (Hanson et al., 2020a, 2025) by assuming that the relative uncertainty stays the same, i.e. the ratio of pre-treatment standard deviation to pre-treatment mean is the same as the ratio of $U_{\mathrm{obs},\mathrm{v},\mathrm{c}}$ to $A_{\mathrm{obs},\mathrm{v},\mathrm{c}}$. The $Q_{\mathrm{obs},\mathrm{v},\mathrm{c}}$ values are the 1σ uncertainty in the linear regression slopes (DeGroot and Schervish, 2018). We chose normalization by uncertainty to ensure the BGNPP_treeshrub variable, which are based on ingrowth core samples that had limited spatial representativeness and only span two years, was not overemphasized compared to the other C fluxes, which have five-year estimates (Hanson et al., 2020a; Iversen et al., 2021b). The normalization additionally ensured that all the variables are unitless, making them intercomparable.

For ELM-MYCI, the total number of preexisting and newly added parameters resulted in an extremely large search space that cannot be covered by 4000 samples. Therefore, we focus on perturbing the most sensitive half (18) of the 37 newly added parameters (Table 1). The most sensitive half was predetermined using one-at-a-time (OAT) perturbation – for each parameter, we generated 50 uniform random samples per parameter between pre-defined upper and lower bounds (Tables S4, 1). The perturbation ranges were based on observations or previous studies when possible (Table S4). The range for q₁₀ of mycorrhizal acquisitions was set wider than the range of q10_hr for HR (Table S3) to explore more values (Table 1). The ranges for the rate constants ($v_{\mathrm{N}/\mathrm{P},\mathrm{froot},j}$, $v_{\mathrm{N}/\mathrm{P},\mathrm{fungi},j}$, $u_{\mathrm{N}/\mathrm{P},\mathrm{fungi},j}$) always span exactly two orders of magnitudes, and were selected based on experimental measurements (Table S5) and similar rate constants in past modelling studies (He et al., 2021; Shao et al., 2023b; Wu and Blodau, 2013). The geometric centre of each range was a manually tuned value that produced acceptable performance at the site during preliminary simulations. The ranges for the half-saturation points ($k_{\mathrm{N}/\mathrm{P},j}$) were set at approximately the modelled annual average soil inorganic NP levels at the SPRUCE site, because experimentally measured values (Table S6) were found to cause PFTs to die in test simulations. The ranges for the C-cost factor of mycorrhizal N acquisition (c_N) were selected based on past observed exchange ratios for EcM (Bogar et al., 2022), assuming ErM has similar behaviour. We did not include the high-end of the observed exchange ratios (>100; Bogar et al., 2022) because those severely restricted mycorrhizal nutrient uptake in the modelling context, possibly due to the 50 % net photosynthesis cap.

We then calculated the standard deviations of $A_{\mathrm{sim},\mathrm{v},\mathrm{c}}$ and $S_{\mathrm{sim},\mathrm{v},\mathrm{c}}$ over all the random samples for each variable and CO₂ treatment. This approach resulted in a 4th order tensor of 37 parameters × 6 variables × 2 CO₂ treatments × 2 coefficients [mean or temperature response slope]. We converted the standard deviation tensor to a same-dimensioned rank tensor by ranking the standard deviations high-to-low across the parameters separately for each variable, CO₂ treatment, and mean or slope. We show the ranks in Fig. S4, with the last three dimensions of the rank tensor compressed into the form of boxplots. The selected most sensitive parameters are those with lowest median ranks in Fig. S4 and are listed in Table S5. We then perturbed those most sensitive 18 parameters in a 4000-member ensemble simulation, calculated the RAE in the same way as done for ELM-OLD, and selected the sample with lowest RAE as the optimized set of parameters.

It is desirable to understand whether our structural change can improve ELM-OLD beyond the capability of parameter optimization. Therefore, we conducted two sets of 4000-member ensemble simulations for ELM-MYCI: one with the preexisting parameters fixed at default levels, and one with the preexisting parameters fixed at optimized levels. Table 2 lists the final four simulations compared in this study. Comparing ELM-OLD_optim to ELM-OLD gives the effect of parameter optimization. Comparing either ELM-MYCI to ELM-OLD, or ELM-MYCI_optim to ELM-OLD_optim, gives the effect of structural modification. Comparing ELM-MYCI_optim to ELM-OLD gives the combined effects of parameter optimization and structural modification.

Table 2 Summary of parameter perturbation experiments and the optimal runs compared in this study.

2.7.2 Verification of parameter constraint

Equation (5) only uses C fluxes but the optimized new parameters are mainly relevant to NP. To examine whether the selected optimal parameter values are indeed from a constrained subregion of the search space by Eq. (5), or simply perform best by chance, we used a clustering metric. This metric, though not using a formal likelihood framework, is consistent with the equifinality concept, where parameters are considered identifiable if the high-likelihood parameter values are concentrated in a small region, and non-identifiable if broadly distributed across the parameter space (Raue et al., 2009). We first normalized all the sampled parameter values to between 0 and 1 using the predefined ranges (Tables S3 and S5). Then, we calculated the Euclidean distances between the normalized parameter values of two types of ensemble member pairs. The first type was all the pairs that can be formed between the best performing 1 % (i.e. 40) members. The second type was all the pairs formed between one value from the best performing 1 % members, and one value from the bottom 99 % (4000 − 40 = 3960) members. If the pairwise distances in the first group are systematically lower than the pairwise distances in the second group, it means the better-performing parameter sets are clustered, which we interpret as being constrained by Eq. (5).

2.7.3 Sobol's sensitivity analysis

While the OAT perturbation helps us to downsize the number of parameters to perturb, it offers limited insight into the full sensitivity of model outputs, because it omits potential interactions among parameters. Using the 4000-member ensembles, we therefore applied a previously established procedure on ELM to calculate Sobol's main-effect and total-effect indices between selected model outputs and all the perturbed parameters (Ricciuto et al., 2018). Sobol's sensitivity indices are based on variance decomposition, where the total variance in a model output is decomposed into fractions attributed to individual parameters and interactions among parameters (Saltelli et al., 2004). The main effect measures the relative fraction attributed to a parameter excluding all interactions. Due to the omission of interactions, the sum of main effects across all the parameters is less than 1. The total effect measures the relative fraction attributed to a parameter including itself and all higher-order interactions involving the parameter. Due to duplicated counting of interactions, the sum of total effects across all the parameters is greater than 1.

3 Results

3.1 Comparison between model and observations

3.1.1 Model performance on C fluxes and leaf area index

On the calibration set of C fluxes, ELM-MYCI_optim best captures the mean and temperature sensitivity of the C fluxes (RAE = 1.54), followed by ELM-OLD_optim (RAE = 1.61), ELM-OLD (RAE = 1.80), and lastly ELM-MYCI (RAE = 1.82; Fig. 2). Also, ELM-MYCI_optim and ELM-OLD_optim have the fewest number of variables falling outside the ± one standard deviation observational uncertainty window (6 out of 28), followed by ELM-MYCI (8 out of 28), and lastly ELM-OLD (9 out of 28; the number of “×'s” of Fig. 2). Among the individual variables, the mean NPP_moss under elevated CO₂, the temperature sensitivity of BGNPP_TreeShrub under ambient and elevated CO₂, and the temperature sensitivity of aggregate NPP under elevated CO₂ are least well-captured, with all four model setups simulating outside the observational uncertainty (Fig. 2a, c, d). ELM-MYCI captures mean HR better than the other three model setups but performs worse on several other variables (mean AGNPP_shrub, temperature sensitivity of AGNPP_spruce under ambient CO₂, mean NPP_moss under ambient CO₂, and temperature sensitivity of NPP_moss under elevated CO₂; Fig. 2). Consistent with our hypothesis (see end of Sect. 1), structural modification improves model performance in capturing the observed large positive response of shrub growth to warming. The temperature sensitivity of AGNPP_shrub under ambient and elevated CO₂ are outside the observational uncertainty in ELM-OLD (Fig. 2c). Using optimized preexisting parameters, ELM-OLD_optim simulates the temperature sensitivity of AGNPP_shrub within the observational uncertainty range, but still underestimates it by ∼ 50 % under elevated CO₂ (Fig. 2c). With structural modification on top of optimized preexisting parameters, ELM-MYCI_optim best captures the temperature sensitivity of AGNPP_shrub, with only 17 %–35 % underestimation (Fig. 2c). Structural modifcation appears to negtatively influence the simulated mean AGNPP_spruce, which is more severely underestimated in ELM-MYCI and ELM-MYCI_optim than in the other two model setups.

Figure 2 Observed versus simulated mean C fluxes ($A_{\mathrm{obs},\mathrm{v},\mathrm{c}}$, $A_{\mathrm{sim},\mathrm{v},\mathrm{c}}$ in Eq. 5) and temperature sensitivity of C fluxes ($S_{\mathrm{obs},\mathrm{v},\mathrm{c}}$, $S_{\mathrm{sim},\mathrm{v},\mathrm{c}}$ in Eq. 5) on the calibration set (2015–2021). RAE is calculated from Eq. (5) for each model setup. The C fluxes names and calculation of the means and temperature sensitivities are explained in Sect. 2.5. The additional “NPP” term was summed up from C fluxes 1–6 in Sect. 2.5. to provide a measure of plot-level primary productivity, but should not be understood as standard NPP because it does not contain coarse root production. The suffix “CO2” means values for elevated CO₂ enclosures and without this suffix means values for ambient CO₂ enclosures. The observational uncertainty intervals are estimated as described in Sect. 2.7.1. The “×” on top of each bar indicates that the simulated value is outside the observed uncertainty interval.

The results of Fig. 2 show aggregated statistics for all the plots. To reveal more detailed causes of model biases, Fig. S7 shows the plot-wise observed and modelled annual time series of the calibration C fluxes. The interannual and inter-plot variabilities are much larger in the observations than all four model setups, especially in 2015, suggesting spatial variability in the real world remains a source of uncertainty in the model-observation comparison (Sect. S1.2). Increase in AGNPP_shrub, decrease in NPP_moss, and increase in HR from the coldest to the warmest treatments are visually discernible in both the observations (Fig. S7c1, d1) and the simulations (Fig. S7c2–c5, d2–d5). Using root mean squared error (RMSE) as a metric on each individual C flux, one can see that ELM-MYCI_optim best captures AGNPP_shrub (Fig. S7c2–c5) as well as the aggregate NPP (Fig. S7f2–f5), ELM-MYCI best captures HR (Fig. S7g2–g5), but the structural modifications did not result in clear improvements on AGNPP_spruce or AGNPP_tamarack (Fig. S7a1–b5). These are consistent with calibration-period performance evaluated on mean and temperature sensitivity (Fig. 2).

On the evaluation set of C fluxes, ELM-MYCI_optim has the best RAE (4.65), followed by ELM-OLD_optim (RAE = 5.03), ELM-MYCI (RAE = 7.22), and ELM-OLD (RAE = 7.31; Fig. 3). RMSE and RAE values calculated on the individual C fluxes show that ELM-MYCI_optim outperforms the other three model setups on AGNPP_spruce, AGNPP_shrub, and NEE_shrubmoss, but not on AGNPP_tamarack. This is consistent with the better performance of ELM-MYCI_optim on shrub than trees during the calibration period (Figs. 2, S7). None of the model setups correctly captures the positive temperature sensitivities of AGNPP_spruce under both CO₂ treatments and of AGNPP_tamarack under the ambient CO₂ treatment in 2022 (Fig. 3a, b). The high RAE on NEE_shrubmoss likely arises from an underestimation of observational uncertainty for this variable, since it was estimated from the RMSE of the gapfilling algorithm (Figs. S4–S5), not considering pre-treatment spatial variability or uncertainty in the chamber- and plot-level biomass measurements used to harmonize model and observations (Sect. 2.6.1). Nonetheless, such bias does not affect the within-variable ranking across the model setups based on RAE or RMSE.

Figure 3 Observed versus simulated annual C fluxes on the evaluation set (2022–2023), showing the magnitude and temperature sensitivity. The C fluxes names are explained in Sect. 2.6 and “gs” in the unit of NEE_shrubmoss denotes growing season. For each subplot, linear regressions relate annual mean air temperature to the C flux, separately for the ambient CO₂ treatment plots and the elevated CO₂ treatment plots. Shaded regions indicate 95 % confidence bands of the regressions (DeGroot and Schervish, 2018). Scatter dots show the original data points only for the observational regressions; the simulated original data points are omitted for readability. Root mean squared errors (RMSE) were calculated between the simulated and observed annual C fluxes within each subplot. RAE were calculated from Eq. (5), for each C flux separately (subplot-level) and for all four compared C fluxes (figure-level). For RAE calculation, the uncertainty of AGNPP_spruce, AGNPP_tamarack, and AGNPP_shrub are as documented in Sect. 2.7.1. The mean uncertainty of NEE_shrubmoss was set equal to the combined RMSE of the gapfilling algorithm on all the observed data points (15.24 gC m⁻² gs⁻¹, combined from Figs. S4–S5). The temperature sensitivity uncertainty of NEE_shrubmoss was equal to the standard deviation of regression slope like the other variables documented in Sect. 2.7.1.

In addition to C fluxes, we used annual maximum LAI, which is not a calibration variable, as a proxy to compare modelled C biomass to observed. Based on RMSE, ELM-MYCI_optim has lower overestimation in spruce LAI and lower underestimation in shrub LAI compared to the other three model setups, although higher overestimation in tamarack LAI (Fig. S8).

3.1.2 Model portrayal of soil nutrients

The observed resin-exchangeable nutrients increase by about two orders of magnitude with warming in all observed combinations of depth, nutrient, and hummock-hollow (Fig. 4). The net mineralization rates increase only by about one order of magnitude in ELM-OLD and ELM-OLD_optim. ELM-MYCI and ELM-MYCI_optim better capture the observed increases (yellow line in Fig. 4a, yellow and red lines in Fig. 4b, e, f, red line in Fig. 4c).

Figure 4 Relationships between normalized plant-available nutrients and annual mean peat temperatures in the observations and models at different depths of hummock and hollow. The plant-available nutrients are represented by annual total resin-exchange nutrients (NH${}_{\mathrm{4}}^{+}$ and PO${}_{\mathrm{4}}^{\mathrm{3}-})$ in observation and annual mean net mineralization rates (NET_NMIN and NET_PMIN) in the model. Annual values across all the years, enclosures, and CO₂ treatments are plotted together. The normalization procedure is reported in Sect. 2.6.4 and reconciles unit difference between observation and model. Lines are least-squares linear regression lines.

The modelled soil NH${}_{\mathrm{4}}^{+}$ and soluble P concentrations, converted to per volume soil water basis to compare with observed porewater values, show mean levels that are ELM-MYCI > ELM-OLD > ELM-MYCI_optim > ELM-OLD_optim (Fig. S9a, c). The ELM-OLD_optim values are closest to observed NH${}_{\mathrm{4}}^{+}$ values but still about one order of magnitude higher. The ELM-MYCI values are closest to observed soluble P values but more than one order of magnitude smaller. The observed NO${}_{\mathrm{3}}^{-}$ mean concentrations are captured more accurately than NH${}_{\mathrm{4}}^{+}$ or soluble P, with ELM-OLD_optim being closest to observations, followed by ELM-MYCI_optim, while ELM-OLD and ELM-MYCI remain within one order of magnitude of the observations (Fig. S9b). These large discrepancies call into question if these observed and simulated quantities are comparable, even after the units are nominally aligned. Like the mean concentrations, the observed and simulated trends differ substantially in magnitudes. ELM-OLD and ELM-MYCI better capture the observed qualitative transition from negative to positive trends NH${}_{\mathrm{4}}^{+}$ towards the warmer enclosures than ELM-MYCI_optim; ELM-OLD_optim did not capture this transition for the ambient CO₂ enclosures (Fig. S9d).

All the model setups overestimate the peat C and N stock in the shallow soil layers (0–30 cm; Fig. 5a, b). In the deeper layers (30–200 cm), ELM-OLD exhibits underestimation that is exacerbated in ELM-MYCI and remedied in ELM-OLD_optim and ELM-MYCI_optim (Fig. 5a, b). All model setups severely underestimate P stock below 30 cm (Fig. 5c). The total simulated soil organic C stock of all the soil layers is about 190 kg C m⁻² by ELM-OLD, 240 kg C m⁻² by ELM-OLD_optim, 140 kg C m⁻² by ELM-MYCI, and 210 kg C m⁻² by ELM-MYCI_optim. This places ELM-OLD_optim outside the observational uncertainty of the estimated 176 ± 40 kg C m⁻² of the S1 bog (McFarlane et al., 2018), probably due to the overestimation in the shallow layers (Fig. 5a).

Figure 5 Observed pre-treatment peat carbon, nitrogen, and phosphorus stocks and the modelled values in the ambient control plot. Error bars on the observations are the ± standard errors reported in the original studies (Griffiths et al., 2017; Salmon et al., 2021). The definitions of acrotelm, mesotelm, and catotelm are in Salmon et al. (2021).

3.2 Behaviour diagnostics of the modified model

3.2.1 Plant nutrient acquisition response to warming and soil inorganic nutrients levels

Simulation results from ELM-MYCI and ELM-MYCI_optim show that the total N acquisition from organic sources by PATH^myc,org across the three vascular PFTs are about 2.5 gN m⁻² yr⁻¹, and the total P acquisition by PATH^myc,org about 0.07 gP m⁻² yr⁻¹, in the unenclosed ambient plot (TAMB; Fig. 6). The PFT-total N acquisition by PATH^myc,org is mainly accounted by spruce (about 50 % of the total spruce N uptake, Fig. 6a) and shrub (30 %–50 % of the total shrub N uptake in TAMB, lower in the warmed enclosures, Fig. 6c). The total organic P acquisition mainly comes from tamarack (>50 % of the total tamarack P uptake; Fig. 6e). Those numbers and percentages are comparable to the coarse pretreatment estimates for the SPRUCE site, which suggests organic N sources account for about 30 % of total plants N acquisition (2.4 gN m⁻² yr⁻¹ out of 7.6 gN m⁻² yr⁻¹), and organic P sources account for a negligible fraction of total plants P acquisition (0.7 gP m⁻² yr⁻¹; Salmon et al., 2021).

Figure 6 Simulated annual mean nitrogen (N) and phosphorus (P) acquisition by the vascular PFTs in each treatment plot during 2015–2021. The acquisition is equal to inorganic nutrient uptake in ELM-OLD and ELM-OLD_optim, and equal to the sum of all three pathways (inorganic nutrient uptake by uncolonized fine roots [PATH^root], inorganic nutrient acquisition via mycorrhizal roots [PATH^myc,inorg], and nutrient acquisition from organic sources via mycorrhizal roots [PATH^myc,org]) in ELM-MYCI and ELM-MYCI_optim.

ELM-MYCI and ELM-MYCI_optim simulates higher total spruce N acquisition through all three pathways, and greater warming-induced increases in total tamarack and shrub NP acquisitions across all three pathways, than ELM-OLD and ELM-OLD_optim (Fig. 6). The low spruce P acquisition simulated by ELM-MYCI_optim, especially under elevated CO₂, can explain the underestimation of AGNPP_spruce by this model setup (Fig. 2a). The large increases in shrub NP acquisition are driven by PATH^root (Fig. 6c, f), consistent with declining fungal colonization rates and rising soil inorganic nutrients under warming (Figs. S9a–c, S10c).

In addition to changing the responses of plants nutrient acquisition to warming (Fig. 6), ELM-MYCI and ELM-MYCI_optim simulate more flexible responses of plants nutrient acquisition to soil inorganic nutrient contents (Fig. S11). In ELM-OLD and ELM-OLD_optim, the uptake always displays a logistic shape where they first increase with soil inorganic N or P and then plateaus. In ELM-MYCI and ELM-MYCI_optim, the total acquisition across all three pathways can show linear relationships (e.g. Fig. S11b, e, f, ELM-MYCI_optim) or logistic shapes that saturate at higher or lower levels (e.g. Fig. S11a, b, ELM-MYCI).

3.2.2 Net ecosystem exchange responses to warming

All model setups simulate a transition from C sink to C source, that is, negative to positive net ecosystem exchange (NEE), with warming (Fig. 7a). Except for the +0.00 °C elevated CO₂ treatment, the C source strength is ELM-OLD > ELM-OLD_optim > ELM-MYCI > ELM-MYCI_optim. NEE is mainly driven by the balance between gross primary productivity (GPP), heterotrophic respiration (HR), and autotrophic respiration (AR), i.e. NEE ≈ − (GPP − AR − HR). All model setups simulate decreasing GPP with warming, especially in the ambient CO₂ plots, while the sum of AR and HR remains relatively stable. The latter's apparent stability is due to cancellation between decreasing AR, consistent with the GPP change, and increasing HR under warming.

Figure 7 Modelled plot-wise mean and temperature sensitivity of the net ecosystem exchange (NEE) and its balance terms: gross primary productivity (GPP), autotrophic respiration (AR), heterotrophic respiration (HR), and the ratios of AR and HR to GPP. The values for AR and HR are stacked, because the two add up to ecosystem respiration, with the darker shade being AR and the ligher shade being HR. The temperature sensitivities are calculated as linear regression slopes between the annual mean NEE, GPP, AR, HR, AR-to-GPP ratio, or HR-to-GPP ratio values against the annual mean air temperatures during 2015–2021 in each plot. The temperature sensitivities have solid bars when they are significantly different from zero at p≤0.05 (two-sided t-test), and otherwise hollow bars.

The low NEE of ELM-MYCI_optim is because it has the lowest fraction of GPP lost to AR, especially in the warmest enclosures, among all the model setups (Fig. 7c). Like ELM-MYCI_optim, ELM-MYCI has a low AR-to-GPP ratio, but its low NEE is likely driven by the low GPP per se (Fig. 7b, c). In contrast, ELM-OLD_optim has lower NEE than ELM-OLD is because of a high GPP (Fig. 7b) and a smaller fraction of GPP lost as HR (Fig. 7d). Structural modification has little effect on HR to GPP ratios (ELM-OLD vs. ELM-MYCI, ELM-OLD_optim vs. ELM-MYCI_optim; Fig. 7d). The temperature sensitivities of NEE, GPP, AR and HR are generally insignificant within individual enclosures, but exhibit consistent increasing, decreasing, or bell-shaped patterns across the enclosure warming levels in all the model setups (Fig. 7e–l).

3.2.3 Nutrient limitations on plant growth

The difference between structural modification and parameter optimization on AR-to-GPP ratios can be better understood by examining the individual components of AR (Fig. 7). In ELM, AR is the sum of excess, maintenance, and growth respiration (AR = XR + MR + GR). XR represents respiration loss due to nutrient limitation, and is calculated as an nonlinearly increasing function of the percentage of plant biomass existing as NSC (Sect. S1.1.2); MR represents respiration for maintaining regular plant metabolism and is approximately proportional to total living biomass; GR is a small and constant fraction of structural growth and therefore of low interest here (Oleson et al., 2013).

At grid level, the lower fraction of GPP lost as AR in the structurally modified model setups than their unmodified counterparts (Fig. 7c) is mainly driven by lower XR-to-GPP ratios, especially under warming (ELM-MYCI_optim vs. ELM-OLD_optim, ELM-MYCI vs. ELM-OLD; Fig. 8a). Interestingly, parameter optimization induces large decreases in the MR-to-GPP ratio that are offset by large increases in the XR-to-GPP ratio (ELM-OLD and ELM-MYCI vs. ELM-OLD_optim and ELM-MYCI_optim), resulting in small net changes (Fig. 8a). This “trade-off” between XR and MR can be explained by their implicit relationship. At a higher XR-to-GPP ratio, the higher nutrient limitation prevents GPP from being assimilated into structural growth, leading to lower biomass-to-GPP ratio and therefore lower MR-to-GPP ratio. By the same logic, lower XR-to-GPP ratio implies a higher biomass-to-GPP ratio, and therefore higher MR-to-GPP ratio.

Figure 8 Partitioning of excess respiration (XR) and maintenance respiration (MR) parts of autotrophic respiration relative to gross primary productivity (GPP) for all PFTs combined or individual PFTs. The ratios are calculated by dividing the 2015–2021 averages of the PFT-total or PFT-specific terms. XR contributions are stacked on top of MR contributions, so the total bar height reflect how much of the GPP flux is offset by total XR + MR.

Compared to the grid-level ratios, structural modification has strong effects on the ratios calculated between PFT-specific XR, MR, and GPP. For spruce, ELM-MYCI and ELM-MYCI_optim have much higher XR-to-GPP ratio and lower MR-to-GPP ratio than ELM-OLD and ELM-OLD_optim, especially in the colder enclosures (Fig. 8b), which indicates a XR-MR trade-off similar to observed at grid level. For tamarack, ELM-MYCI_optim has much lower XR-to-GPP ratio than ELM-OLD_optim in all the enclosures (Fig. 8c). For shrub, ELM-MYCI and ELM-MYCI_optim have much lower XR-to-GPP ratio than ELM-OLD and ELM-OLD_optim in the warmest enclosures (Fig. 8d). As such, spruce and shrub are likely the main drivers behind the more rapid decline of the grid-level XR-to-GPP ratio in the structurally modified models than ELM-OLD or ELM-OLD_optim (Fig. 8a). The weaker XR-MR trade-off in tamarack and shrub likely reflects the lower importance of MR in the two PFTs compared to spruce.

3.2.4 Nutrient limitations on heterotrophic respiration

Nutrient limitation of HR occurs in the default and modified models when the available soil inorganic N or P cannot satisfy the total demand of the plants and immobilization demand from soil litter and organic matter decomposition (Eq. S4, Sect. S1.1.3). The immobilization demand arises mainly in the decomposition step from plant litter to SOM (Oleson et al., 2013; Schimel and Bennett, 2004). The C : N and C : P in the plant litter pool depend on litter inputs and are usually higher than the C : N and C : P of the downstream SOM pools, which are fixed parameters (Oleson et al., 2013). As a result, the process immobilizes additional soil inorganic NP to meet the C : N and C : P of the SOM pools. Because ELM-MYCI and ELM-MYCI_optim allow plants to access the NP in plant litter pools via mycorrhizal roots (Sect. 2.2), the C : N and C : P the litter pools increase, resulting in greater immobilization demand per unit decomposition, which is proportional to HR. Consistent with this expectation, for each unit of HR, the corresponding actual immobilized inorganic N in ELM-MYCI and ELM-MYCI_optim are higher than in ELM-OLD_optim and ELM-OLD at the same level of HR (Fig. S12a). The same effect is not seen in P (Fig. S12b). Instead, ELM-OLD_optim exhibits considerably higher P-immobilization per unit HR than the other three model setups. All the model setups exhibit sensitivity to warming, which suggests P-limitation on HR is more affected by litter quality changes created by relative changes in primary productivity among the PFTs than N-limitation on HR (Fig. 2).

3.3 Parameter sensitivity analysis

3.3.1 Constraint of model parameters

The distance metric (Sect. 2.7.2) shows the top-performing 1 % parameter values are statistically significantly closer to each other (smaller distances) than to the remaining 99 % parameter values (larger distances) in all three ensemble simulations (Table 2, Fig. 9). The significant separation means the C fluxes can constrain the preexisting and newly added parameters. The distances are least well-separated for ELM-MYCI_ENS (Fig. 9b), which uses the same un-optimized parameters as ELM-OLD for the unchanged model processes. Those suboptimal parameter values may have caused biases that the new model processes cannot compensate for, leading to unstable optimized values in the new parameters in the ELM-MYCI_ENS simulations.

Figure 9 Probability density distributions of the Euclidean distances between the parameter values of two groups of pairs of ensemble members: between the best-performing 1 % members (shaded bars), and between the best-performing 1 % members and the other 99 % members (solid line). Displayed probability densities are pooled from all pairs in each group. The p-values in each panel indicate whether the two groups are significantly different, using two-sided t-test for the mean values of two independent samples.

3.3.2 Sensitivity of model outputs to parameter values

We show the parameter sensitivities of ELM-MYCI_optim_ENS in the main text (Fig. 10), since the finding of Sect. 3.3.1 suggests this may be a more reliable perturbed parameter ensemble than ELM-MYCI_ENS (Fig. S11), and the sensitivities of preexisting parameters (Fig. S12) are similar to the well-reported findings of past studies (Meng et al., 2021; Ricciuto et al., 2018).

The relative sensitivity of model outputs to each parameter is approximately the same whether assessed using total effects or main effects (compare the rows in Fig. 10). The grid total GPP is sensitive to the parameters of all three vascular plants and grid-level parameters (topmost bar in each panel of Fig. 10). The grid total NEE and vegetation C (TOTVEGC) are most sensitive to spruce parameters, especially the sensitivity of fungal colonization rate to soil inorganic N (b_j, Eq. S12) and maximum rate of inorganic N acquisition via mycorrhizal association ($v_{\mathrm{N},\mathrm{myc},j}$, Eq. S23). The grid total HR and total soil organic C (TOTSOMC) are most sensitive to shrub parameters, especially the maximum rate of N uptake via fine root ($v_{\mathrm{N},\mathrm{froot},j}$, Eq. S33). The C variables in any vascular PFT are mainly determined by the parameters specific to that PFT, especially the maximum uptake/acquisition rates of NP ($u_{\mathrm{N},\mathrm{myc},j}$ and $u_{\mathrm{P},\mathrm{myc},j}$ for spruce (Eq. S17), $v_{\mathrm{P},\mathrm{froot},j}$ for tamarack (Eq. S33), and $v_{\mathrm{N},\mathrm{froot},j}$ for shrub (Eq. S33)). Moss C variables are sensitive to the parameters of all three vascular plants and grid-level parameters.

Figure 10 Sobol's main-effect and total-effect sensitivity indices of selected C-balance variables to the newly added model parameters, calculated from ELM-MYCI_optim_ENS. For better display, the indices are partitioned into subpanels according to whether it is a PFT-specific or column-level parameter. Stacking the bars across the four panels in each row gives the sum of the main or total effects over all the perturbed parameters, which are also displayed as a grey line for reference in each panel. The C-balance variables in each panel are grouped according to whether it is a column-level, spruce, tamarack, shrub, or moss variable. Parameter definitions can be found in Table 1 and equations referred therein.

Compared to ELM-MYCI_optim_ENS, the sensitivities derived from ELM-MYCI_ENS are similar for tamarack and shrub parameters (Figs. 10b, c, f, g, S11b, c, f, g). For spruce parameters, the model outputs of ELM-MYCI_ENS are more sensitive to the fungal colonization rate to soil inorganic N (b_j, Eq. S12 and Fig. S11a, e), whereas ELM-MYCI_optim_ENS are more sensitive to the maximum organic NP acquisition rates via mycorrhizal association ($u_{\mathrm{N},\mathrm{myc},j}$ and $u_{\mathrm{P},\mathrm{myc},j}$ for spruce, (Eq. S17; Fig. 10a, e). Still, the two sets of parameters fulfil similar functions, with the former controlling all the nutrient uptake/acquisition pathways (see M_myc,j in Eqs. S17, S23, and S33) and the latter only controlling the organic pathway (Eq. S17). For column level parameters, ELM-MYCI_ENS is mainly sensitive to the Q₁₀ of NP acquisition rates (q₁₀, Eq. S13) and has little sensitivity to the C cost of mycorrhizal nutrients acquisition to the plant (c_N, Eqs. S25, S30) (Fig. S11d, h). ELM-MYCI_optim_ENS exhibits the same contrast, albeit less strongly (Fig. 10d, h).

Comparing between the newly added (Figs. 10, S11) and preexisting parameters (Fig. S12), one can see the newly added parameters exhibit more inter-PFT interactions. That is, the C variables of each PFT are even more strongly determined by the parameters specific to that PFT in ELM-OLD_optim_ENS (Fig. S12a, b, c, e, f, g) than in ELM-MYCI_optim_ENS (Fig. 10a, b, c, e, f, g) or ELM-MYCI_ENS (Fig. S11a, b, c, e, f, g). Additionally, the grid-level C variables are mainly responsive to spruce parameters in ELM-OLD_optim_ENS (Fig. S12a, e), compared to both spruce and shrub parameters in ELM-MYCI_optim_ENS (Fig. 10a, c, e, g) and ELM-MYCI_ENS (Fig. S11a, c, e, g). Moss C variables are more responsive to spruce parameters in ELM-OLD_optim_ENS (Fig. S12a, e), compared to shrub parameters in ELM-MYCI_optim_ENS (Fig. 10c, g) and ELM-MYCI_ENS (Fig. S11c, g). The main and total effects of ELM-OLD_optim_ENS are close to 1 when summed over all the parameters (Fig. S12 grey lines), while the main effects of ELM-MYCI_ENS and ELM-MYCI_optim-ENS are much smaller than 1 and the total effect much greater than 1 (Figs. 10, S11, grey lines). The larger difference between total and main effects in the modified models means the newly added parameters have stronger statistical interactions than the preexisting parameters in the default model.

4 Discussion

4.1 Summary of evaluation performance and remaining gaps

We present a development of the ELMv2-SPRUCE model to replace the photosynthesis-driven, inorganic-only plant nutrient uptake with three pathways that consider influences from fine root biomass, fungi-colonization level, and plant access to organic nutrients through fungi-colonized roots (Sect. S1.1). Although EcM and ErM are only implicitly represented and simplifications are made in the treatment of the C cost of mycorrhizal acquisition and the organic NP sources (Sects. 2.2, S1.1.8 and S1.1.9), the modified model shows improved performance. Compared to parameter optimization only (ELM-OLD_optim), structural modification (ELM-MYCI_optim) leads to lower RAE and RMSE on vegetation C fluxes, especially in the shrub-moss community (Figs. 2, 3), improved LAI in two out of the three modified PFTs (Fig. S8), improved qualitative similarity to resin-exchange nutrients (Fig. 4), and similar performance on soil total C-N-P stocks (Fig. 5). Interestingly, structural modification imposed on unoptimized preexisting parameters (ELM-MYCI) does not improve RAE (Fig. 2) and the new parameters appear ill-constrained (Fig. 9b). These findings indicate the biases in ELM-OLD arise from both incorrect parameter values and inadequate process representation (Bastrikov et al., 2018). The strong performance of ELM-MYCI_optim, after optimizing newly-added parameters on top of pre-optimized preexisting parameters, supports stepwise calibration as a viable strategy for land surface models when the parameter search space is large (Ma et al., 2024).

One notable finding is that ELM-MYCI_optim captures the observed large positive warming response of AGNPP_shrub better than the parameter-optimized default model (ELM-OLD_optim; Figs. 2b, 3c). This increasing growth is accompanied by rising NP uptake via PATH^root and unchanging NP acquisition from organic sources via PATH^myc,org (Fig. 6c), consistent with our initial hypothesis and previous finding at the Mer Bleue peatland site (Shao et al., 2023b) that declining dependence on ErM drives shrub growth under warming. Recent analysis of minirhizotron data at the SPRUCE site shows increasing specific root length for the shrubs with deeper water tables as a consequence of warming (Weber et al., 2025a, b). This shift towards more acquisitive fine root trait (Bergmann et al., 2020; Weber et al., 2025a, b) is not yet considered in the current study and might partially explain the remaining underestimation in the temperature sensitivity of AGNPP_shrub. The simulated high importance of direct fine root uptake in shrubs at the SPRUCE site differs from the simulated > 90 % dependence on fungal-mined organic N for shrubs at Mer Bleue (Shao et al., 2023b). This wide range is comparable to past observations (Hilman et al., 2024; Hobbie and Hobbie, 2006; Yin et al., 2022) and might reflect inter-site difference, wherein SPRUCE is a more southern site with lower shrub fractional cover and higher porewater inorganic N compared to Mer Bleue (Kennedy et al., 2018; Shao et al., 2023b).

ELM-MYCI_optim more severely underestimates mean AGNPP_spruce than the other model setups (Fig. 2a), likely because it simulates stronger P limitation on spruce (Fig. 6d). The stronger P limitation, in turn, may be because the modelled peat P stock and soil inorganic P levels are generally too low (Figs. 5c, S7c). With the enhanced shrub growth in ELM-MYCI_optim, the remaining inorganic P becomes insufficient to support spruce growth (Fig. 2a). All model setups underestimate mean AGNPP_spruce in the ambient CO₂ plots, fail to capture the reversal from negative to positive response to warming over time (Fig. 3a), and fail to capture the observed lack of response to elevated CO₂ (Fig. 2a). Those biases might be related to seasonal variations and acclimation in spruce photosynthetic parameters (Dusenge et al., 2024; Jensen et al., 2019), gradual acclimation in tree respiration or hydraulics (Hanson et al., 2025), delayed response to elevated CO₂ or increased C allocation belowground (e.g. to roots, mycorrhizal fungi, or exudates; Duchesneau et al., 2024a; Norby et al., 2010, 2024; Palmroth et al., 2006). Addressing them will require future process developments.

The severe underestimation of the temperature sensitivity of BGNPP_treeshrub in all model setups (Fig. 2b) may be due to high uncertainty in the ingrowth core observations or the presence of dynamic above-to-belowground allocation in response to warming (Drewniak, 2019; Rehschuh et al., 2022). Dynamic allocation is not yet in the model processes of ELM-OLD or ELM-MYCI. Note the ingrowth core observations used in Fig. 2 only span 2016–2017, while the other evaluation variables span 2016–2019 and 2021 (Hanson et al., 2020a). The uncertainty problem will be remedied as additional years of ingrowth core and minirhizotron observations are completed for the SPRUCE site. The initial ingrowth core observations do suggest the fine root biomass of the trees and shrub are more sensitive to warming than their aboveground NPP (compare the fine root biomass reported in Fig. S1 of Malhotra et al. (2020a) to the observed aboveground NPP values this paper Fig. S7).

The persistent bias in NPP_moss and its temperature sensitivity is unsurprising. ELMv2-SPRUCE cannot yet simulate the decline of moss growth with warming (Norby et al., 2019), instead depending on assigned observed fractional covers (Sect. 2.4). Missing processes may include shading from shrub, inaccurate photosynthesis-water relations, and microbial relationships (Carrell et al., 2019; Norby et al., 2019; Petro et al., 2023; Shi et al., 2021). The persistent bias in vertical distribution of peat C and N may be due to insufficient vertical mixing of the soil decomposition pools (Oleson et al., 2013). The persistent bias in HR and peat P stock might be due to inaccurate C : P in the SOM pools (Figs. 2b, 5), and the current lack of consideration of fungal respiration (Sect. S1.1.7). Because the HR-to-GPP ratio is remarkably invariant to our current structural modification (Fig. 7), focused parameter investigation and structural modification on the soil decomposition model may be needed to address the HR bias.

The discrepancy between observed and modelled porewater concentrations (Fig. S9) suggests the values are not directly comparable, which may be due to missing process representation of the adsorption of NH${}_{\mathrm{4}}^{+}$ to inorganic and organic matter surfaces (Eick et al., 1999; Matschonat and Matzner, 1996), inaccurate partitioning between labile P and soluble P (Yang et al., 2023), and underestimation of peat P stock (Fig. 5). Better matches between the temperature sensitivities of normalized model NP mineralization and normalized resin-available NH${}_{\mathrm{4}}^{+}$ and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (Fig. 4) suggest the model is better at capturing relative changes in plant nutrient availability than absolute sizes.

4.2 Impact on ecosystem productivity

We found that parameter optimization reduces the strength of NEE increase under warming (Fig. 7a) via higher GPP and lower HR (Fig. 7b). The GPP effect is likely explained by the dominant control of the photosynthetic parameter “flnr” (fraction of leaf N in in Rubisco enzyme) on grid- and PFT-level C balances in ELM-OLD_optim_ENS, and the HR effect likely by the parameter “q10_hr” (Q₁₀ for heterotrophic respiration; Fig. S14). Interestingly, the overall NEE balance is most strongly affected by the Q₁₀ parameter of spruce MR (Fig. S14), despite this parameter having little effect on the other C-balance terms (Fig. S14) – this might be a case of emergent phenomena (Brient, 2020; Wang et al., 2022) and worth future modelling and empirical investigations.

Structural modification reduces the extent of NEE increases under warming (Fig. 7a) via lower AR (Fig. 7c), which is driven by decreases in XR with warming (Fig. 8a), especially in spruce and shrub (Fig. 8b, d). The XR-MR trade-off at grid-level and in spruce (Sect. 3.2.3) mitigates the XR-driven decreases in AR, demonstrating a case of nonlinear feedback. The greater declines in XR in ELM-MYCI and ELM-MYCI_optim compared to ELM-OLD and ELM-OLD_optim are directly due to declines in NSC (Sect. S1.1.2), which is in turn likely due to reduced nutrient limitation (Fig. 6a, c, d, f) under warming. Although NSC is additionally affected by the C cost of the mycorrhizal pathways (Sect. S1.1.8), the C cost should decrease with warming as the importance of the mycorrhizal pathways decline (Fig. 6a, c, d, f). Therefore, this mechanism cannot explain the simulated NSC decline with warming. The lower nutrient limitation under warming implies greater increase of plant carbon use efficiency (CUE), especially for the shrub PFT where most of AR is due to XR (Fig. 8d). This modelling result is consistent with empirical evidence that EcM-tree association is key to explaining the negative correlation between CUE and latitude in northern boreal regions (Mäkelä et al., 2022). This consistency supports the effectiveness of our implicit approach as a parameter-efficient framework. Biologically, the CUE-latitude correlation is driven by complex interactions between plant C flow to mycorrhizal colonization level, and nutrient availability (Mäkelä et al., 2022; Shao et al., 2023b).

In the structurally modified model setups, SOM decomposition becomes more N-limited because of the acquisition of N from plant litter pools via PATH^myc,org (Fig. S12, Sect. S1.1.9). Surprisingly, the higher N-immobilization per unit HR only corresponds to slightly lower mean HR in ELM-MYCI and ELM-MYCI_optim (Figs. 2b, d, 7d, h), meaning the lower availability of organic N in plant litter pools are offset by higher inorganic N uptake. This pattern suggests that decomposition suppression through nutrient competition – classically termed the Gadgil effect when mediated by EcM fungi (Fernandez and Kennedy, 2016) – may be limited in strength under our simulated conditions. Although our model does not explicitly represent fungal guild interactions, the modest reduction in decomposition is consistent with studies showing that such effects are highly context-dependent and often confounded by litter quality, niche partitioning among EcM, ErM, and saprotrophs, and priming processes (Fanin et al., 2022; Mielke, 2022; Shao et al., 2023a). P-limitation on immobilization appears to be controlled by more complex factors than N-limitation (Fig. S12b) and may be affected by current model bias in peat P stock (Fig. 5). Additionally, soil decomposition process in ELMv2-SPRUCE do not explicitly simulate microbial biomass and guilds, and therefore might misrepresent the partitioning of immobilization demand between the external nutrient uptake and the cycled nutrient between dead and live microbial biomass, which is especially important for P (Duchesneau et al., 2024a; Schmidt et al., 1997). Given this limitation in the ELMv2-SPRUCE soil decomposition model and current model limitation in treating organic NP sources (Sect. S1.1.9), the HR results should be interpreted with caution.

4.3 Limitations, Uncertainties, and Future directions

4.3.1 Limitations of the ELM-MYCI equations

In ELM-MYCI, the partitioning of plant nutrient demand between non-mycorrhizal and mycorrhizal pathways is a linear function of soil inorganic N concentration (Eq. 1). The focus on N is likely acceptable for northern boreal regions, where plant growth is N limited or co-limited by N and P, but extending the model into temperate or tropical regions will require considering P limitation (Du et al., 2020). A second question is whether the covariate, soil inorganic nutrient concentration (Eq. 1), is simplistic and should be replaced with plant nutrient status. Plant exudates and phytohormones play a large role in the early establishment of EcM and AM, making plant nutrient status a mechanistically defensible choice (Garcia et al., 2015), but this mechanism is much less established in ErM (Genre et al., 2020). EcM colonization may also respond to soil inorganic N directly, because assimilating inorganic N is more C-expensive than assimilating organic N (Garcia et al., 2015). Therefore, the relative appropriateness of the two covariates will be best determined via multi-site, multi-PFT simulations across N and P gradients. The current ELM-MYCI equations already have two terms related to plant nutrient status: the N-limitation and NSC availability multipliers (Eqs. 3–4). Therefore, to implicitly encode an effect on colonization rate, one can simply let the parameter value in the N-limitation multiplier (Eq. S15) differ between PATH^root and the two mycorrhizal pathways.

The current formulation for PATH^myc,org makes a few significant simplifications that do not conform to real world processes (Sect. 2.3). The restriction of mycorrhizal organic NP acquisition to plant litter pools may underestimate total mycorrhizal organic NP acquisition. This, in turn, could lead to underestimation of mineral NP limitation on HR and overestimation of the NP limitation on plant growth, contributing to the biases shown in Fig. 2. The acquisition rate does not vary with soil organic NP content or the recalcitrance of the accessed plant litter pools (Eq. 4), and the upper bound (0.0001 ⋅ total plant litter pool size per hourly time step) is more a sanity constraint rather than a realistic upper bound on fungi-accessible organic matter. This may result in underestimation of the changes in acquisition rates across environmental gradients. In reality, mycorrhizal nutrient acquisition involves two sequential steps: depolymerization followed by uptake of small organic molecules for N (Talbot and Treseder, 2010), and enzymatic hydrolysis and uptake of released phosphate for P (Plassard et al., 2011). Each step is describable by Michaelis-Menten kinetics (Näsholm et al., 2009; Plassard et al., 2011). In the organic matter-rich peatland environment, the first step would be much more limited by fungal enzyme availability than substrate availability, making a colonization-proportional formulation (Eq. 4) acceptable. Extending the model to mineral soil environments will require revisiting the formulation. Although the rate of the first step may be higher for labile than recalcitrant organic matter (Talbot and Treseder, 2010), current understanding of mycorrhizal enzyme activity is largely from a biochemical perspective, i.e. characterizing individual compounds, rather than from an ecological perspective, i.e. classifying numerous compounds into measurable fractions that significantly differ in rates (Lavallee et al., 2020; Oleson et al., 2013). Therefore, the current even distribution of organic NP acquisition across the plant litter pools, proportional only to pool sizes, serves as a starting point that can be refined as ecological-scale data become available. Modelling the second step in Michaelis-Menten form would require explicit representation of the turnover of mono- and oligomeric fractions of organic N pools, or of soluble P pools in microsites, which are beyond the scope of this study.

The total mycorrhizal acquisition of inorganic and organic NP is capped by 50 % of net photosynthesis per time step (Sect. 2.3). In the real world, belowground C allocation to mycorrhizal fungi peaks later in the growing season than aboveground photosynthesis (Högberg et al., 2010). Therefore, the current model formulation may underestimate mycorrhizal NP acquisition later in the growing season. In principle, this could lead to underestimated NP limitation on the heterotrophic microbes later in the growing season and overestimated NP limitation during leaf expansion in the following spring. The practical effect in the current ELMv2-SPRUCE, however, may be minor. The model also does not yet separate heterotrophic microbes and mycorrhizal fungi in soil decomposition. The model also uses fixed C : N : P ratios for plant structural tissues, and always evaluates nutrient limitation based on concurrent soil nutrient availability rather than nonstructural nutrient storage pool sizes.

4.3.2 Uncertainty of the ELM-MYCI parameters and potential constraints

The new equations in ELM-MYCI introduces a fairly large number of new parameters (in total 37 in Tables S4 and 1) compared to ELM-OLD. Although the annual C fluxes demonstrated ability to constrain these parameters (Fig. 9), the ability may have benefited from the temperature gradient, with implied inorganic nutrient gradient, across the SPRUCE treatments. Extending ELM-MYCI to other sites will therefore require more direct constraints.

The highest priority should be given to uptake rate constants ($v_{\mathrm{N}/\mathrm{P},\mathrm{froot},j}$, $v_{\mathrm{N}/\mathrm{P},\mathrm{fungi},j}$, $u_{\mathrm{N}/\mathrm{P},\mathrm{fungi},j}$) and half-saturation points ($k_{\mathrm{N}/\mathrm{P},j}$) because of the high sensitivity of model outputs to these parameters (Sect. 3.3.2). Many experimental measurements of these parameter exist, but the data points are concentrated in temperate regions and crop species, and exhibit high cross-species uncertainty (Craig et al., 2025; Griffiths and York, 2020; Tables S5–S6). Therefore, these data points are likely useful for large-scale simulations, but less so for site-scale simulations involving only a few species or simulations targeting underrepresented ecosystems like peatlands. In those cases, total plant NP uptakes, estimated from paired NPP and whole-plant C : N and C : P measurements, provide a reasonable nutrient-centred constraint. The efficiency of the parameter search can be increased in future studies by constraining the search space with well-validated qualitative findings in experimental data – for example, mycorrhizal roots generally have higher rate constants for inorganic NP on a per unit fine root biomass basis than uncolonized roots (Craig et al., 2025; Plassard and Dell, 2010).

Other large sources of uncertainty in the modified model include the fungal colonization fractions (Fig. S10), relative contributions of different pathways to total plant NP acquisition (Fig. 8), and the transfer of plant C to mycorrhizal fungi (Sect. S1.1.8). EcM colonization is known to respond to large temperature, water, and nutrient gradients and covary with fine root traits (Ostonen et al., 2011, 2017; Xie et al., 2021). Stable isotope fractionation studies (δ¹⁵N, δ¹³C) provide the most direct constraint on mycorrhizal contribution to plant demand and plant C allocation to mycorrhizal fungi (Hawkins et al., 2023; Hilman et al., 2024; Hobbie and Hobbie, 2006; Yin et al., 2022). Observations of all these quantities vary between 0 %–100 % (Hawkins et al., 2023; Hilman et al., 2024; Hobbie and Hobbie, 2006; Ostonen et al., 2011, 2017; Xie et al., 2021; Yin et al., 2022) and root tip counts at SPRUCE corroborate the high uncertainty (Fig. S10). Therefore, multi-site model-data integration will be needed to capture the broad pattern and prevent overfitting.

Root morphology parameters come from measurements in this study (root diameter – r_j and root tissue density – ρ_j in Table S5). These two traits belong to, respectively, the collaboration and conservation gradients in the root economics space, and both are correlated with other traits (Bergmann et al., 2020; Craig et al., 2025). Thinner roots tend to depend less on fungal partners and have higher uptake rate constants. Lower root density tends to be associated with higher root N content, turnover rates, uptake rate constants, yet higher half-saturation points. These qualitative understandings, along with global root trait observations from, e.g. the Fine-Root Ecology Databases (Iversen et al., 2021a), provide good constraints for ELM-MYCI in future multi-site or regional extensions.

4.3.3 Other potential future work

Several process developments in ELM are being carried out during the same time as this work, including division of fine root biomass into transport, absorption, and mycorrhizal pools (Wang et al., 2023), nutrient-responsive dynamic above-belowground allocation (Knox et al., 2024), and explicit microbial biomass pools and dissolved organic matter dynamics (Ricciuto et al., 2021). These developments can be merged with this work to relax many of the assumptions listed in Sect. 4.3.1 and improve the model results. Explicit simulation of fungal and heterotrophic microbial biomass will enable separating mycorrhizal fungal from non-mycorrhizal respiration and the modelling of fungal necromass, which has different decomposability from saprotrophic residues due to higher melanin content, particularly in ErM fungi (Fernandez et al., 2019). Other efforts of interest may include separating the behaviours of EcM and ErM, spruce photosynthesis and MR, Sphagnum growth, and considering direct plant uptake of small organic N molecules (Näsholm et al., 2009). Nonetheless, process-development should consider the limited constraint available from empirical observations (Fig. 9) and ensure the complexity is commensurate with our ability to check model accuracy and interpret cause and effect in model responses.

The modified ELMv2-SPRUCE has extended capabilities compared to the default model, e.g. fine tune nutrients competition relationships between PFTs using the maximum uptake/acquisition rates and half-saturation parameters (Tables S4–S5), assimilating nutrient uptake kinetics data, and testing the ecosystem impacts of changing fine root traits. Those improvements will enable new hypothesis testing and more accurate modelling of peatland C, N, and P cycling. It will also be interesting to compare the modified ELMv2-SPRUCE with other models that use similar fine root-based uptake rules (Knox et al., 2024; Zhu et al., 2019) and/or have mycorrhizal representations (He et al., 2018; Shao et al., 2023b; Sulman et al., 2019), and to test the performance of the model at multi-site to regional scale and its implications for carbon cycle feedbacks to the climate system.

5 Conclusions

We present a development on the ELMv2-SPRUCE model to replace default, photosynthesis-driven nutrient uptake processes with fine root and implicit mycorrhizal pathways, allowing more realistic processes like the access to organic nutrients by mycorrhizal roots and the dependence of plant nutrient uptake on fine root biomass, fungi colonization level, and environmental conditions. The modified ELMv2-SPRUCE model better captures the observed large increase in shrub growth under whole ecosystem warming than the default model, as well as overall measured C fluxes and resin-exchange nutrients response to warming. The modelled increase in shrub growth is accompanied by stable fungi-mediated nutrient acquisition from organic matter, and several fold increase in direct fine root inorganic uptake, supporting our initial hypothesis that the observed increase in shrub growth is driven by a shift from mycorrhizal outsourcing to direct fine root uptake strategy. Non-validated comparisons to the default model show the modified model simulates less nutrient limitation on plant growth under warming, resulting in weaker C-sink to C-source transition, and more flexible relationships between plant nutrient acquisition and soil inorganic nutrient concentrations. Outstanding model biases and caveats indicate needs to improve non-mycorrhizal processes for spruce and Sphagnum moss growth, and above-to-belowground allocation. Other future developments may add or refine representation fine root trait responses to warming, more realistic organic nutrient access, shifts in allocation, and mycorrhizal fungal biomass turnover. Overall, the new model is a useful tool for model-data integration, hypothesis testing in ecosystem carbon-nitrogen-phosphorus cycling, and investigating boreal peatland responses to environmental change.