Evaluation of simulated N cycling using observations from a ¹⁵N tracer experiment in a mixed deciduous forest

Nitrogen availability constrains terrestrial carbon uptake and storage, yet large uncertainties remain in the magnitude of the effect, because the interactions of the carbon and nitrogen (N) dynamics are challenging to observe in undisturbed ecosystems at relevant timescales. Long-term experiments with ¹⁵N tracer applications allow the study of the nitrogen cycle in a fairly undisturbed manner, and they are therefore a valuable data source to test the biogeochemical dynamics simulated by terrestrial biosphere models. In this study we applied the model QUINCY (QUantifying Interactions between Terrestrial Nutrient CYcles and the climate system), which includes an explicit representation of terrestrial ¹⁵N fluxes and pools. We used observations from a long-term (10-year) ¹⁵N tracer experiment in a temperate deciduous forest to evaluate the nitrogen dynamics simulated by QUINCY. Recovery in soil N dominated overall ecosystem ¹⁵N recovery in both observations and simulations over the long-term. The observed gradual movement of the ¹⁵N tracer to lower soil layers was also captured by the model. However, in the short-term the modeled uptake and losses of ¹⁵N for leaves and fine roots were too fast, and recovery in litter and surface soil was too slow, indicating that the model likely overestimated plant competitiveness for newly added N relative to soil microbes. Downward vertical transport of ¹⁵N tracer in the soil was slower in the model compared to measurements, which may be indicative either of too little bioturbation or vertical transport via leaching. Overall, the QUINCY model results showed good agreement with the observations, making it a valuable tool for studying long-term nitrogen dynamics. Running the model over an extended period indicated that the ecosystem retained a very large share of the added ¹⁵N tracer (>90 %), and that this retention persisted over multi-decadal timescales. This study shows that explicit inclusion of isotopic tracers allows for a more thorough evaluation of carbon-nitrogen turnover and dynamics and thereby can contribute to reduce uncertainties in modelling nitrogen cycling and constraints in terrestrial ecosystems.

1 Introduction

Climate change influences plant-atmosphere interactions as ecosystems experience new conditions (Canadell et al., 2022). It is crucial to have a better understanding of the capacity of vegetation to take up atmospheric CO₂ and the fate of terrestrial carbon (C) storage, as potential feedbacks may attenuate or further accelerate climate change. The bioavailability of nitrogen (N) plays a pivotal role in the carbon sink potential of vegetation and soils, as bioavailable N is scarce in many ecosystems owing to limited inputs from atmospheric deposition and biological nitrogen fixation (LeBauer and Treseder, 2008), but plants require N for the the photosynthetic machinery and for constructing new biomass (Bonan, 2016). Understanding the role of the N cycle in the terrestrial C cycle is especially important as terrestrial N limitation is expected to intensify in response to the rising atmospheric CO₂ (Cambron et al., 2025; Rogers et al., 2017).

At the same time, human activities have not only perturbed the C cycle, but also the natural N cycle through industrial combustion processes and agricultural fertilization (Gruber and Galloway, 2008). As a result, the amount of reactive N has increased leading to increased N deposition across many regions of the world (Galloway et al., 2008). These anthropogenic influences on the nitrogen cycle have diverse effects on natural ecosystems, such as eutrophication, air pollution and biodiversity loss (Gong et al., 2024; Schlesinger, 2009). Increased N deposition also influences C sequestration, by stimulating tree growth and altering soil decomposition processes (Zaehle and Dalmonech, 2011). The magnitude of the effect of N deposition in forests has been debated (Gurmesa et al., 2022; Magnani et al., 2007; Sutton et al., 2008; Reay et al., 2008). The global N-induced increases in forest net carbon uptake have been estimated to range from 130 to 345 Tg C yr⁻¹ by several studies using different approaches (de Vries et al., 2014; Du and de Vries, 2018; Fleischer et al., 2019; Jain et al., 2009; Thomas et al., 2010; Schulte-Uebbing and de Vries, 2018; Wang et al., 2017; Zaehle and Dalmonech, 2011), with a smaller recent estimate of 41 Tg C yr⁻¹ from a data-based approach (Schulte-Uebbing et al., 2022).

There has been a trend in declining N availability in unfertilized ecosystems globally, which has lasted for decades (Craine et al., 2018; Mason et al., 2022; Bassett et al., 2026). Declining N deposition resulting from air pollution legislation has contributed to reduced N inputs to forests in Europe and north-eastern North America (Lloret and Valiela, 2016; Schmitz et al., 2019). Declining N deposition, taking place simultaneously with decreases in sulphur deposition, has been shown to have positive effects on some forest ecosystems (Dalton et al., 2024).

An important uncertainty in the response of the carbon cycle to N additions is the fate of the added N in the ecosystem (Nadelhoffer et al., 1999b). The response depends on the fraction of N that is retained in the ecosystem rather than lost to leaching or gaseous emissions, and the partitioning of the retained nitrogen among pools with a narrow C:N ratio, such as soil organic matter, or a wider C:N ratio, such as woody biomass (Gurmesa et al., 2022; Nadelhoffer et al., 1999b). Quantifying the fate of the incoming N to the ecosystem is therefore crucial.

Terrestrial biosphere models (TBMs) are useful tools for predicting the trajectories of terrestrial C stocks in soil and vegetation (Friedlingstein et al., 2025; Seiler et al., 2022), because they can, in principle, account for the numerous interactions between the N and C cycles. However, modelling of the N cycle is challenging because many processes are involved but comprehensive data sets to evaluate the models are scarce (Meyerholt et al., 2020), partly because many aspects of the N cycle are difficult to measure (Vicca et al., 2018). Results from manipulation experiments, including both N fertilization and tracer addition studies, are important resources for the further evaluation and development of the N cycle description in TBMs (Thomas et al., 2015).

N fertilization experiments in forests have a long history, the earliest dating back to the early 20th century (Högberg et al., 2017). These experiments have been valuable for estimating optimal N fertilization for forest growth as well as for increasing our understanding of the N cycle in these ecosystems (Högberg et al., 2017). These data have also been used for model development (Meyerholt and Zaehle, 2015; Caldararu et al., 2020). However, these experiments have limitations, as fertilized ecosystems are no longer in their natural state, and fertilization may lead to N saturation or species replacement (LeBauer and Treseder, 2008; Van Houtven et al., 2019). In model development it is often preferable to have observations of ecosystems at their natural state. Additionally, it is difficult to trace the effects of individual N cycle processes beyond their influence on tree biomass, although such insights are needed to further improve our understanding these processes and their representation in models.

Isotopic tracer experiments can overcome these issues, as very small additions of nitrogen in the form of highly enriched ¹⁵N are sufficient to allow the tracer to be tracked across different soil and vegetation pools (Buchmann et al., 1995; Lin et al., 2024; Nadelhoffer et al., 1999a; Schleppi and Wessel, 2021). Traditionally, these observations are challenging and expensive to obtain. However, improved measurement techniques have made such observations more common, providing more data (Schleppi and Wessel, 2021) that can also be used for model development. These observations have revealed many important aspects of the nitrogen cycle that would otherwise have been difficult to study (Craine et al., 2015; Robinson, 2001). For example the uptake of N by the forest canopy has been quantified (Ferraretto et al., 2022; Li et al., 2025).

These isotopic tracer experiments have been used to test models without an explicit ¹⁵N cycle (Cheng et al., 2019; Thomas et al., 2013) and have provided insights into nitrogen turnover rates and competition for N between plants and soil microbes. A more direct use of these ¹⁵N observations is to compare them with the simulated ¹⁵N cycle in models that explicitly represent it, as in the TBM QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) (Thum et al., 2019). While site-level models used for studying tracer experiments have been developed (e.g., Inselsbacher et al., 2013; Müller et al., 2014; Liu et al., 2026) and some simple models do describe the ¹⁵N cycle (Vitousek et al., 2024), QUINCY, as a TBM has the advantage of having fully coupled energy, water, carbon, and nutrient cycles, allowing the study of several aspects of biogeochemical cycling.

QUINCY has already been used by Caldararu et al. (2022) to explain the drivers of long-term trends in natural-abundance measurements of ¹⁵N in plant tissues. In the present study, we compared QUINCY with observations from a ¹⁵N tracer experiment conducted in a temperate mixed deciduous forest, where the fate of the applied ¹⁵N tracer has been tracked for a decade. We specifically aimed to evaluate the temporal dynamics of the tracer signal in different pools over an 11-year period, the partitioning of the ¹⁵N signal among different ecosystem compartments, how well the model reproduced these patterns, and what this assessment reveals about potential improvements to nitrogen cycle modelling. We also conducted a simulation experiment to investigate the fate of the ¹⁵N over decadal timescales and to consider the impact of N inputs on the carbon pools of the ecosystem. In particular, we address the following research questions:

• How realistic are the simulated turnover rates of plant N pools?

• How well does the model simulate retention of the ¹⁵N tracer in the ecosystem?

• How long does the ¹⁵N tracer remain in the ecosystem and where does it ultimately end up?

2 Materials and methods

2.1 ¹⁵N experiment and the experiment site

2.1.1 Site description

The Arnot Forest (42°17^′ N, 76°38^′ W, 503 m elevation) is located in Central New York State, USA. The mean annual temperature is 7.8 °C and the annual mean precipitation is 930 mm (Goodale, 2017). The forest consists mostly of second-growth mixed hardwoods established naturally after harvest in between 1873 and 1887 and fires in 1900–1911 (Fain et al., 1994). The dominant species include white ash (Fraxinus americana), red maple (Acer rubrum), quaking aspen (Populus tremuloides), basswood (Tilia americana) and with smaller amounts of eastern hemlock (Tsuga canadensis) (Goodale, 2017). Tree biomass is approximately evenly divided between ectomycorrhizal and arbuscular mycorrhizal species. The soil consists of 18 % clay, 41 % silt and 40 % sand and the pH is 4.4, measured at the 0–10 cm depth (Fahey et al., 2013). Leaves usually emerge in early to mid-May and most leaf fall occurs by the end of October.

2.1.2 Tracer addition experiment

In 2007, 0.21 kg ha⁻¹ of ¹⁵N tracer was applied as 99 at. % ¹⁵N−KNO₃, split into three equal applications of 0.07 kg ha⁻¹ applied on 30 April, 31 July and 30 October (Goodale et al., 2015). The total amount of the three applications was less than 3 % of one year of nitrogen deposition at the site. The tracer was sprayed onto the forest floor. Measurements of ¹⁵N in fine and coarse roots from the uppermost soil layer 10 cm of soil were conducted in late November 2007. Tree N and ¹⁵N pools (foliage, wood, bark, coarse and fine roots) and soil pools down to 50 cm were measured in summers 2006 and 2008, one year before and after the tracer additions. Additional measurements were made in 2012–2013 (5–6 years after tracer addition), when soil and root samples were taken in October 2012 and the aboveground pools were sampled in June 2013. The final measurements of the plant and soil ¹⁵N pools were collected in 2017, 10 years after the tracer addition. The wood pools were not measured at this time. In the observations, up to 26 samples were taken for the foliage pool, and around 20 for other aboveground plant pools and soil pools (Goodale et al., 2015; Goodale, 2017). Fine and coarse roots were represented by up to 10 samples, with somewhat smaller sample sizes in deeper layers. The standard error for the stock estimates were around 5 %–10 %, with highest uncertainties for the fine roots. The ¹⁵N observations were on the order of 10 %–20 %, again with the highest standard error for the fine roots.

The fraction (%) of the added ¹⁵N tracer recovered in each of the pools ${{}^{\mathrm{15}}\mathrm{N}}_{\text{Rec}}$ was computed from each pool's N storage (M_pool N) and the change in mean atom ¹⁵N following tracer addition (${{}^{\mathrm{15}}\mathrm{N}}_{\text{post}}$) relative to the mean pre-addition value (${{}^{\mathrm{15}}\mathrm{N}}_{\text{pre}}$) and the total mass of added tracer ($M_{{}^{\mathrm{15}}\mathrm{N}\text{added}}$):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{{}^{\mathrm{15}}\mathrm{N}}_{\text{Rec}}& \left(\mathit{\%}\right)=\\ & {\displaystyle \frac{(\text{at.\hspace{0.17em}\%}{{}^{\mathrm{15}}\mathrm{N}}_{\text{post}}-\text{at.\hspace{0.17em}\%}{{}^{\mathrm{15}}\mathrm{N}}_{\text{pre}})\cdot M_{\text{pool\hspace{0.17em}N}}}{M_{{}^{\mathrm{15}}\mathrm{N}\text{added}}}}\cdot \mathrm{100}\end{array}\end{array}$$

2.2 Model description of QUINCY

QUINCY is a terrestrial biosphere model that simulates the energy, water, carbon, and nutrient (nitrogen and phosphorus) cycles (Thum et al., 2019). For simplicity, in this study we ran the model with soil phosphorus availability set to non-limiting levels, such that plant and soil biogeochemical processes operated at average observed P contents. Alongside the nitrogen cycle, an explicit representation of the ¹⁵N cycle has been implemented, largely following the description of isotopic discrimination processes by Robinson (2001). Different ecosystems are represented by different plant functional types (PFTs) and in the site-level version each site is represented by a single PFT.

We outline some of the basic characteristics of the model here, and a full description of the model is provided by Thum et al. (2019). Plants comprise three short-lived pools (leaves, fine roots and fruits), a non-structural storage pool (reserve) and a seasonal, non-respiring pool (labile), along with three structural tissue types (sapwood, heartwood and coarse roots). All vegetation and soil pools in QUINCY contain carbon, nitrogen and ¹⁵N. Photosynthesis in the multilayered canopy is calculated following the formulation of Kull and Kruijt (1998). Photosynthetic parameters depend on the leaf N concentration (Friend et al., 1997). Leaf biomass and N concentration respond to soil N availability (Hyvönen et al., 2008). Plant nutrient uptake depends on fine-root biomass, soil mineral nutrient concentration and plant nutrient demand.

Newly acquired carbon, nitrogen and ¹⁵N are first allocated to the labile pool, from which they can be used for tissue growth, respiration or storage in the reserve pool. The potential growth rate is influenced by the PFT-specific allometric relationships among leaves, sapwood and fine and coarse roots. Belowground carbon allocation increases under nutrient or water stress. Actual growth is further controlled by the nutrient demand of each plant pool and the nutrients available for growth, and may also be limited by air temperature and soil moisture. Flexible stoichiometry is permitted for leaves and fine roots.

Symbiotic biological nitrogen fixation (BNF) by microbial–tree partnerships is described as a dynamic trade-off between carbon and nitrogen opportunity costs, following Rastetter et al. (2001) and Meyerholt et al. (2016). Asymbiotic BNF in the litter layer depends on soil temperature, soil moisture, and the N deficit in litter decomposition. Nitrogen fixed through this mechanism is added to the fast soil organic matter pool. The parameters governing BNF were adjusted to align with the findings of Davies-Barnard and Friedlingstein (2020).

The simulated soil has a vertically layered structure, with 15 layers of exponentially increasing soil thickness down to 9.5 m. The soil biogeochemical processes largely follow the CENTURY model approach (Parton et al., 1993). Each soil layer contains organic pools consisting of metabolic (soluble), structural (polymeric), and woody litter pools, fast- and slow-overturning organic matter (SOM), as well as inorganic ammonium (NH₄) and nitrate (NO₃) pools. The litter and SOM pool turnover rates depend on soil temperature and moisture, described by first-order kinetics. The stoichiometry of the litter pools is determined by the respective plant pool stoichiometries along with the plant resorption capacity for each nutrient. The fraction of resorption to the labile pool for leaves before shedding is 50 % and for wood 20 % during conversion of sapwood to heartwood. Fine-root mortality is assumed to be dominated by predation and therefore no nutrient resorption occurs from fine roots. The C:N stoichiometry of the fast SOM pool is affected by available inorganic N while the slow SOM pool stoichiometry is assumed to be fixed. Plants and soil microorganisms compete for the inorganic nutrients according to their nutrient demand. Nitrification and denitrification processes depend on the aerobic status, temperature, and moisture of each soil layer (Zaehle and Dalmonech, 2011). As a further development of the baseline version of the model (Thum et al., 2019), nitrate assimilation by soil microorganisms (i.e. immobilisation of nitrate from the solute pool to the fast SOM pool) was implemented. Although this allows more rapid uptake of ¹⁵N from nitrate in the soil, this change has limited overall impact on the model results, as N immobilisation by microbial growth is primarily determined by the stoichiometric imbalance between litter and SOM, as well as by the potential litter decomposition rate in the model.

Isotopic discrimination of ¹⁵N occurs during biological nitrogen fixation, ammonification, plant and microbial N uptake, and processes associated with nitrification and denitrification. Fractionation parameters were taken from Robinson (2001) (Table A1). However, the model does not represent fractionation associated with mycorrhizal-mediated N uptake. Biomass growth and N retranslocation are assumed not to fractionate N.

2.3 Model set-up

In this study, the model was run at the site scale, with the site represented by a summergreen broadleaf deciduous forest PFT. QUINCY uses a half-hourly time step and requires air temperature, humidity, precipitation, air pressure, short- and longwave radiation, wind speed, atmospheric CO₂, ¹³CO₂, ¹⁴CO₂ mole fractions and NH_x, NO_y and PO₄ deposition rates as forcing variables. Additional site information includes geographical coordinates and soil physical and chemical properties (texture, bulk density, rooting depth and soil depth, inorganic P content). The turnover rates of vegetation and soil pools for this PFT are listed in Table S1 in the Supplement.

We simulated the Arnot Forest using meteorological data from the CRU JRA v2.1 dataset (University of East Anglia Climatic Research Unit and Harris, 2020), which were disaggregated to half-hourly resolution using a statistical weather generator (Zaehle and Friend, 2010). Annual atmospheric CO₂ concentrations were taken from Le Quéré et al. (2018) and N deposition data from Lamarque et al. (2010, 2011). The prescribed N deposition (0.87 $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. S1 in the Supplement) was similar to the observations at the site (0.7–0.9 $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). The level of N deposition decreased by about one-third between 2000 and 2020, from 1.3 to 0.9 $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at a long-term monitoring site operated by the US Environmental Protection Agency, which is located close to the site (Butler et al., 2015; U. S. Environmental Protection Agency Clean Air Markets Division, 2024). This level of N deposition is considered to represent a transition from high to moderate deposition (Xie et al., 2024). Soil texture and pH were taken from the observations (Fahey et al., 2013). We did not include woody litter in the analysis, because it plays only a minor role in the δ¹⁵N cycle over the timescale considered here. We combined observations of new wood, old wood and bark and compared them with simulated sapwood and heartwood, labile, and reserve pools. Although bark can be considered an active pool in the ¹⁵N cycle (Goodale, 2017), it is not explicitly represented in the model.

First, a 1000-year spinup was performed by cycling driving data between the years 1901 and 1930. After this a model simulation was performed for 1901–2018, using transient climate, CO₂ concentrations and N deposition for the site. In the simulations the ¹⁵N tracer was applied as NO₃ solute in the uppermost soil layer on the same dates as the tracer application was performed in the experiment. For the scenario run extending further 30 years beoynd 2018, the climatic forcing was generated using the present-day climate and the CO₂ from the RCP8.5 scenario (Riahi et al., 2011).

Additional sensitivity tests were performed to assess the extent to which the simulation results depend on the assumed fraction of N retention during leaf senescence, as well as on the turnover rate of soil organic matter pools. In these tests, we varied the coefficient determining the amount of nutrients returned to the reserve pool before leaf senescence, or increased the turnover rates of the soil pools, thereby reducing their size. As these changes affected the overall model performance without improving agreement with the observed characteristics of the carbon and nitrogen cycle or strongly affecting the temporal evolution of the tracer distribution, we used the default parameters for this PFT rather than performing rigorous parameter tuning.

3 Results

3.1 Comparison between magnitudes of observed and simulated pools and fluxes

First, we investigated the overall model performance in reproducing the observed characteristics of the Arnot site – vegetation and soil pools as well as C and N fluxes – because differences in pools between the observations and simulations influence the δ¹⁵N values from the tracer experiment (Table 1). If a pool was overestimated in the simulations, the magnitude of the ¹⁵N signal would be biased low given observation-based flux partitioning. It is important to note that the model was not calibrated to this particular site and the default parameter values were used.

Table 1 Observed (in bold) and simulated carbon and nitrogen pools and fluxes. The simulated values are 20 year averages. Soil value refers to the top 50 cm of soil.

The simulated annual gaseous N fluxes were of the same order of magnitude as the observations, particularly considering the large uncertainty range in the observations. Although key indicators of productivity were comparable to the observations (wood production and leaf C, which is related to leaf area index), the model showed notable deviations from the observed partitioning of carbon and nitrogen in vegetation. The model underestimated the amount of belowground carbon, especially fine-root biomass, while generally exhibiting a tighter C:N stoichiometry in leaves, fine roots, and woody biomass. Discrepancies in leaf litterfall were possibly the result of under-reporting in the observations combined with the omission of the coniferous hemlock species in the model. Comparison of leaf mass and leaf litterfall stoichiometry suggests that leaf N resorption during leaf senescence was lower than the assumed 50 % in the model. In the observations, the mortality exceeded the wood production, which was not reproduced by the model. Instead, the model estimated that wood production and mortality were roughly balanced. The N flux associated with mortality was larger than the observed, probably because the observation-based estimate of mortality only accounted for the mortality of woody biomass, whereas the mortality term in the model also accounted for leaf, root, and fruit biomass.

The litter (also denoted as Oi) C pool was overestimated in the simulations compared with the observations. Given the discrepancy between observed litter C:N ratio and litterfall C:N ratio, it appeared likely that processes during litter decomposition already affected observed litter C:N, whereas the model's litter pool represents fresh, undecomposed litter. The simulations strongly overestimated soil C, as it was almost twice the observed value (Table 1), where the fast-turnover pool contributed to 17 % of the total soil C pool and 18 % of the total soil N pool. The soil C:N ratio of the combined simulated fast and slow soil pools was similar to the observations. The observed soil C:N ratios declined strongly with soil depth, whereas, as assumed by the model, the simulated soil C:N ratio remained nearly constant across soil layers (Table S2 in the Supplement).

In the observations (Table S3 in the Supplement), the natural abundance δ¹⁵N values of leaves (−2.1 ‰) and coarse roots (−2.4 ‰) were most depleted and the wood (−1.2 ‰) was the most enriched vegetation pool. The simulated leaf δ¹⁵N values (−2.3 ‰ vs. −2.1 ‰) and wood (−1.2 ‰ vs. −1.1 ‰) were similar to the observations. The simulated δ¹⁵N in coarse roots (−1.4 ‰) was less depleted than in the observations (−2.4 ‰). The natural abundance δ¹⁵N levels in the litter layer and soil were more enriched than in the simulations, with the difference being largest in the deepest soil layers (Table S3).

3.2 δ¹⁵N signal: Temporal evolution in vegetation

In the experiment the tracer was sprayed on top of the litter layer. In the simulations, we added the incoming tracer to the N deposition flux, from which it directly entered the soil-solution NO₃ pool of the uppermost soil layer. This slight difference may cause some discrepancies in the comparison, as in the model, the soil-solution NO₃ pool of the first soil layer is directly accessible to plants and soil microbes. Given the timescale of this analysis, however, we assumed that delays in tracer transfer from the litter layer into the soil would not substantially affect tracer accessibility to plants and soil microbial biomass, and would therefore have little effect on our results.

Similar to observations, simulations showed a rapid δ¹⁵N enrichment of both leaves and fine roots in the year after the tracer addition, but the simulated decline of δ¹⁵N enrichment over the next decade was faster and larger than observed for these pools (Fig. 1). Simulated leaf δ¹⁵N enrichment (Fig. 1a) showed a small increase already in 2007, but reached its highest level of 51.4 ‰ in 2008 (hereafter, all values reported are relative to the initial natural abundance δ¹⁵N), the year after tracer addition. One year later, the simulated δ¹⁵N enrichment decreased by more than 50 % and reached a low, steady enrichment (5.9 ‰) by 2012. The observations showed a smaller enrichment peak than the simulations (27.1 ‰), but the subsequent decline was slower and less pronounced than the simulated decline.

Figure 1 The temporal change of the observed (in black) and simulated (in red) δ¹⁵N from the reference value in the different vegetation pools (a – leaf, b – fine roots, c – coarse roots, d – wood, including sapwood, heartwood and bark) and the sizes of the corresponding nitrogen pools (e). Note: The line representing leaf δ¹⁵N is not continuous due to the absence of leaves during the dormant season.

The fine root (Fig. 1b) observations showed a similar behavior to that of leaves, with a large increase in δ¹⁵N enrichment (46.6 ‰) during the year of labeling and some decrease in the first year after labeling (27.1 ‰), followed by a smooth decline by approximately 50 % over the following decade. The first observed δ¹⁵N value in 2007 includes only the uppermost layer (0–10 cm), whereas the other observations include fine roots in all the soil layers. The value in the upmost layer is larger than averaged across all soil layers, as it is where the largest amount of fine root biomass is located and the δ¹⁵N enrichment decreases with depth. The simulations approximately matched the observed δ¹⁵N enrichment peak (with a value of 25.5 ‰) observed in 2008, then showed a faster and larger decline than in the observations, dropping by two-thirds within five years and stabilizing at a lower level of δ¹⁵N enrichment compared to the observations.

The observed coarse-root δ¹⁵N values (Fig. 1c) showed a smaller peak in 2007 (33.1 ‰) than fine roots, but had a larger δ¹⁵N enrichment peak in 2008 than either the leaf or fine root pools. This suggests that the tracer was first captured by fine roots and then was subsequently transferred to the coarse roots in the following year in the observations. After peaking in 2008, the observed values declined, levelling off at approximately half of the δ¹⁵N peak by 2012. The simulations did not show such a distinct peak, only a small, gradual increase in the first four years (up to 6.6 ‰) after which the simulated δ¹⁵N remained nearly constant, decreasing only slightly to 6.1 ‰, well below the observed value (16.1 ‰).

The observed wood δ¹⁵N (including bark) increased by 4.9 ‰ in 2008 compared to the pre-labeling value and increased slightly further by 2012 (Fig. 1d). The simulations showed a larger δ¹⁵N enrichment peak (8.1 ‰), which declined to 2.7 ‰ within four years.

3.3 δ¹⁵N signal: Temporal evolution in soil

The observed soil δ¹⁵N enrichment after 2006 was much larger than the simulated one in all soil pools (Fig. 2), which may partly result from the overestimation of the soil N pools in the simulation in both shallow (3.4-fold overestimated in simulation) and deep soil (2.2-fold overestimated in the simulation) (Fig. 3, Table S2). The observations showed a δ¹⁵N enrichment peak of 106.4 ‰ in the litter layer compared with the pre-experiment value and a decrease to about 10.4 ‰ over the following decade. The simulated ¹⁵N enrichment peaked later and only reached a value of 28.9 ‰ and a declined to a level of 4.6 ‰ over the following decade. In the observed shallow soil layer (0–10 cm), the tracer caused a δ¹⁵N enrichment of 9.2 ‰ in 2008 and a slow increase to a level of 12.4 ‰ in a decade. By contrast, the simulations showed an increase only of 3.2 ‰ and the δ¹⁵N enrichment started to slowly decrease, reaching 2.7 ‰ after a decade. The deeper soil layers experienced an increase in δ¹⁵N by 1.7 ‰ in 2008 and this value slowly increased up to 2.1 ‰ in a decade. The simulation, however, showed an increase of only 0.43 ‰ in a decade.

Figure 2 The temporal evolution of the observed (in black) and simulated (in red) δ¹⁵N signal from the reference value before labeling in the different soil pools (a – litter, c – soil between 0–10 cm depth, e – soil between 10–50 cm) and the sizes of the corresponding measured and simulated nitrogen pools (b, d, f). The dashed line shows in (a) the soluble litter only, in (b) the fast soil organic matter pool between 0–10 cm depth and (c) the fast soil organic matter pool between 10–50 cm.

When only the soluble litter was considered (Fig. 2a) the agreement between simulation (simulated peak: 40.9 ‰) and observation was slightly better. The reason for this discrepancy lies in the immediate “adsorption” of the tracer signals when sprayed onto the litter surface, which cannot be represented by the model, since uptake of inorganic N by the litter pool is not a flux represented in the model. Instead, the δ¹⁵N tracer was added directly to the soluble NO₃ pool, from which it was available for plant and microbial uptake, the latter of which generates fast SOM, as well as leaching and gaseous emissions. In the simulation, the δ¹⁵N signal in the litter originated solely from litterfall of enriched plant material, and thus most of the added δ¹⁵N tracer was immobilised by plants or soil microbes. In fact, while Fig. 2b and c show that the model did not capture the magnitude of the bulk SOM δ¹⁵N signal, when only the fast SOM pool was considered, the observed change in SOM was well captured in the topsoil layer (peak enrichment of 12.2 ‰) (Fig. 3b). The observed δ¹⁵N values in the deeper layers showed little variation, similar to the simulated changes, although a slight propagation of the surface signal to lower layers can be seen in later years.

Figure 3 The temporal evolution of the observed δ¹⁵N at different depths (a) and simulated (b) and the soil profiles of nitrogen (organic and inorganic) in observations and simulations (c).

The vertical distribution of the ¹⁵N tracer showed that the observed soil was more enriched in 2006 than the simulated soil throughout the depth profile (Fig. 3). Furthermore, the observed δ¹⁵N signal in the soil was most depleted in the top layer, and most enriched at around 0.25 m depth. In the simulations the soil became steadily more enriched towards deeper layers in 2006.

The temporal changes in soil δ¹⁵N were much larger in the observations than in the simulation, as already shown in Fig. 2. In 2008 the increase in observed soil δ¹⁵N extended to a depth of 30 cm, and during the following years the enrichment continued to increase slowly in the upper soil layers. In the simulation, the temporal dynamics were different. The top layer was most enriched in 2008 and subsequently became depleted. In 2008, the signal had not yet reached the layers below 20 cm in the simulation. During the following years, the enrichment resulting from the δ¹⁵N tracer addition gradually propagated to deeper soil layers. Comparison of the observed and simulated soil N profiles showed that the simulation overestimated the soil N pool in the top layer, whereas below this layer the vertical gradient was more stable and similar to that observed (Fig. 3).

3.4 Recovery of tracer in present day climate and future

The total recovery of the tracer was close to 100 % in the simulation (Fig. 4), whereas it was around 70 % in the observations. Both the observations and the simulation showed a large fraction of the δ¹⁵N signal remaining in the shallow soil layer, and the relative contribution of deeper soil layers increased only slowly over time. The observations showed a noticeable fraction of the δ¹⁵N signal remaining in the litter layer, which was almost depleted by 2017, whereas in the simulation the litter pool was of minor importance.

Figure 4 The temporal evolution of the observed (a) and simulated (b) recovery rates of different plant and soil pools in different years.

The δ¹⁵N signal in coarse roots was pronounced in the observations, but not in the simulation. Both the observations and the simulation showed a contribution from the fine roots. The simulation showed a substantial contribution from the wood pool already in 2008, which remained substantial thereafter. In the observations the transfer of δ¹⁵N to the wood pool occurred in 2008 and increased in 2012 (wood was not measured in 2017; values for 2012 were retained for visualization). This difference occurred partly because we attributed all labile and reserve ¹⁵N to the woody pool, whereas in reality the signal may be partitioned between coarse roots and woody reserve storage. The simulation showed a large fraction of the recovered δ¹⁵N in leaves in 2008 (9.5 %), while the observations showed a recovery of only 1.9 % for that year. When soil and plant pools were considered separately, the observations showed a prominent signal in the soil already in 2008, and its relative contribution decreased over the following years (Table 2). Recovery of the label in the plant pools peaked at 11 % the year after the tracer addition, then stabilized at 8 % in years 5 and 10. In the simulations, the relative contribution to plants was larger than in the observations in 2008, but was similar to the observations in 2012 and 2017 (Table 2).

Table 2 The recovery rates in the soil and vegetation pools in different years. The observed values are in bold and the simulation results in a regular font. Also the percentual amount of the whole recovered signal in the plants is shown.

In order to track how much of the δ¹⁵N signal remained in the ecosystem over longer timescales, we made a scenario simulation for an additional 30 years. The simulated total recovery rate of the system remained high until the end of the simulation, suggesting a relatively closed simulated nitrogen cycle at the site (Fig. 4b). Most of the signal remained in the soil, with the deeper soil layers gaining δ¹⁵N signal from the top soil layers. Over time, the woody pool became an increasingly important vegetation pool for storing the δ¹⁵N signal.

The vertical profile of simulated δ¹⁵N in soil (Fig. 3b) became more uniform over time, while the uppermost soil layers became depleted. A comparison of the simulated δ¹⁵N values in different ecosystem pools in 2017 and 2048 (Table 3) indicates that all vegetation pools, except the wood pool, showed a gradual decrease in δ¹⁵N during this period. The δ¹⁵N in the litter was slowly decreasing. In the soil, δ¹⁵N decreased in the fast SOM pool but increased in the slow SOM pool.

Table 3 Simulated δ¹⁵N recovery in different plant and soil pools in 2017 and 2048.

4 Discussion

4.1 Temporal dynamics of the ¹⁵N tracer signal in vegetation

Our simulation showed how the ¹⁵N tracer signal applied to the forest floor was transferred among ecosystem pools over decadal timescales. Despite some conceptual differences in the description of the manipulation experiment, the QUINCY model successfully simulated many features observed in the observations. The comparison with the observations was especially interesting when evaluating the N turnover rates of the vegetation pools. These suggest that the model's turnover rates of N were too fast for the leaves and fine roots. This occurred despite the seemingly higher N retention of leaves during senescence in the model. It is important to note that the model underestimated the C and N pool size of fine roots and C pool size of the leaves, and this might have some influence on the simulated dynamics. Although the simulated total plant C pool (12 771 g C m⁻²) was close to the observed value (11 140 g C m⁻²), the other vegetation C pools, except for wood, were underestimated. This suggests that QUINCY has some difficulties in representing C allocation. The simulated total N in vegetation (86.3 g N m⁻²) was over double the observed amount (39.2 g N m⁻²), reflecting the overestimation of soil N.

Since this is a deciduous forest and the turnover time of fine roots is set to 0.7 years in the model, it is likely that the ¹⁵N tracer was transferred too quickly away from the compartments that are used to build leaves and fine roots. In the observations, the coarse roots showed a large ¹⁵N peak, whereas the simulation showed only a gradual increase. The current model structure does not accommodate a quick response from the coarse roots, that currently have a turnover time of eight years. The observations showed a larger signal in the coarse roots, because the movement of the signal to the upper tree parts required transport, e.g., via diffusion or xylem transport, and additionally the coarse-root observations included xylem tissue (Tegeder and Masclaux-Daubresse, 2018; Van Der Heijden et al., 2015). In QUINCY, we could account for the larger xylem concentration in coarse roots by partitioning the labile and reserve pools among coarse roots, leaves and woody biomass, whereas in this study for simplicity we have assumed they were allocated to sapwood. Given the simulated peak ¹⁵N signal in labile+reserve pool, partitioning 11 % of this pool to coarse roots would reduce this mismatch. This suggests that the timescale of the transfer of the ¹⁵N signal within the plant (from coarse roots to woody and leaf compartments), is potentially significantly larger (on the order of months to few years), than the nearly instantaneous redistribution of nitrogen in the plant as assumed by the model.

Finally, it should be mentioned that the definition of coarse and fine roots differed slightly between the model and the observations. In the observations, roots thicker than 1 mm were considered coarse roots (Goodale, 2017), whereas the threshold was 2 mm in the model. The 1 mm cut-off point for the fine roots in these observations was chosen to align with earlier studies conducted at this site (Goodale et al., 2015; Goodale, 2017). It has also been shown that fine roots with diameters less than 1 mm constitute the majority of fine-root biomass and exhibit the greatest sensitivity to nutrient availability in this and similar forest types (Fisk et al., 2004; Frey et al., 2025). We therefore do not expect this discrepancy between the simulation and observations to have a significant impact on the comparison.

In the model, the tracer was added to the top layer of solute NO₃ which enabled direct access by plants and bypassed the litter layer. The trees took up this nitrogen into labile and reserve plant N pools in the model. This caused the simulated plants to have more of the ¹⁵N than the soil than observed at the end of 2007. However, by 2008, the observations showed a relatively larger fraction of the tracer in the plants than the simulations. This difference suggests that the simulated plant pools had too high turnover rates for nitrogen, as was also seen separately for the leaves and fine-root compartments. In addition, the strong uptake of the tracer by the litter layer and followed by the gradual release of ¹⁵N could also have influenced the plant signal in the observations. In the observations, the recovery of the ¹⁵N in coarse roots remained substantial during the last years of the experiment and this supports the hypothesis that the tracer was still being slowly released from the litter and being transported to the plants. Generally the relative recovery of ¹⁵N between plants and soil was comparable between the simulation and the observations in later years. However, during the few first years the model estimated nearly twice the recovery of ¹⁵N in plants (∼21 %) as the observations (∼11 %), in line with what has been found out in the earlier modeling studies that modeled general N dynamics rather than explicit ¹⁵N cycling (e.g. Cheng et al., 2019).

4.2 Temporal dynamics of the ¹⁵N tracer signal in soil

Implementation of nitrate assimilation by soil microorganisms was a key process for successfully simulating the tracer experiment at the Arnot Forest. The observed and measured ecosystem N pool magnitudes differed, which also influences the interpretation of ¹⁵N temporal dynamics. These differences were pronounced in the soil, with the modeled total soil N pool (1656 g N m⁻²) was more than twice as large as the observed pool (744 g N m⁻²). Since in the model the slow soil pool was large, but it was not taking up a lot of the ¹⁵N signal due to its lower turnover rate, its large size caused low values of the ¹⁵N tracer signal. This effect was demonstrated by isolating the fast pools (Fig. 2), that better captured the observed changes in the magnitude of the ¹⁵N signal in the uppermost soil organic matter layer.

In the observations, the ¹⁵N signal reached deeper soil layers, down to 0.35 m (observed layer depth 0.3–0.4 m), after only one year and at 0.45 m (observed layer depth 0.4–0.5 m) within five years. In the simulation, the downward movement of the ¹⁵N signal was slower. Several factors may have contributed to this finding, including a too large uptake of the signal by vegetation and soil microbes, a too small movement of NO₃, or an underestimation of the transport by bioturbation. Bioturbation in the model is represented as a simple diffusive flux as in Koven et al. (2013), but with the rate declining with soil depth driven by a decline on the fine root biomass density. This decline might be too strong in QUINCY, but also the amount of C allocation to roots was underestimated in the model, which might have contributed to a too weak bioturbation flux. A further contribution may be that a fraction of the leaching flux was simulated to occur via lateral flow and not deep drainage, effectively partitioning the leaching loss into a lateral component that prevented a fast and deep penetration of the ¹⁵N signal in the soil.

In the observations, the top soil layers became enriched as the ¹⁵N from the overlying litter layer was incorporated into the soil organic matter. The top layers showed some depletion by 2018. The simulation showed larger changes in the top layers over time and also the ¹⁵N tracer was also transported to deeper soil layers, which resulted in a flattening of the vertical ¹⁵N soil profile (Fig. 3). These differences between the model results and observations stem from conceptual differences. It is likely that in the observations the ¹⁵N tracer got attached to the litter layer and was being released slowly still after several years, resulting in a continued enrichment of the top soil layer. This is possible in the real world because the observed litter pool includes early stages of decomposition, typically associated with an increase in N concentration due to immobilisation of mineral N by decomposers, which facilitates the rapid incorporation of ¹⁵N. In the model, the litter pool only contains fresh litter and any derivates, including immobilized mineral N, that are transferred to the fast SOM pool. Disentangling these conceptual differences from any model structural deficiencies is challenging. It is of note that the considerable overestimation of the soil N pool that occurred across all soil depths partly contributed to the simulated soil ¹⁵N signal being more depleted than observed.

The ¹⁵N tracer was applied as NO₃ remained for a very long time in the soil according to both the observations and the model analysis. This is in line with other natural forest ecosystem ¹⁵N studies (Lin et al., 2024; Templer et al., 2012), where the largest compartment of recovered ¹⁵N was also the soil organic matter. It has been suggested that the NO₃ can be incorporated into soil organic matter through abiotic processes (Fitzhugh et al., 2003). These kinds of processes are not yet implemented in the current model set-up.

4.3 N status of the forest

The natural-abundance soil δ¹⁵N differed between the observations and the simulation in 2006, prior to the start of the experiment, with observed soil δ¹⁵N values being more enriched than those in the simulation. Greater soil δ¹⁵N enrichment is often associated with increased N losses via gaseous pathways which strongly discriminate against δ¹⁵N (Craine et al., 2018; Feng et al., 2023). This could indicate a less closed N cycle and weaker N limitation in the observations than in the simulation. However, following the initial tracer loss, the observed δ¹⁵N recovery remained relatively stable over time, suggesting relatively efficient long-term N retention and a relatively closed N cycle. The mycorrhizal fungi present at the site are known to strongly fractionate N isotopes (Hobbie and Ouimette, 2009; Hobbie and Högberg, 2012), which may have also contributed to the enriched soil δ¹⁵N values in the observations. Leaf stoichiometry can be considered a proxy for leaf N status (Sheng et al., 2021). The simulation showed a lower leaf C:N ratio than the observations, suggesting weaker N limitation in the simulation.

N deposition declined during the study period. In N-limited ecosystems such a decline might lead to a more closed N cycle. However, in boreal forests, the role of increasing atmospheric CO₂ has been shown to be a more important driver for a reduction in N availability than changes in N deposition (Bassett et al., 2026) and a similar response may occur in temperate forests. Overall, changes in N limitation during the 11-year study period are unlikely to have influenced our conclusions.

4.4 Recovery rates and future

The ecosystem-scale recovery of the ¹⁵N tracer is very challenging (Craine et al., 2015). Recovery of the tracer remained close to 100 % in the model throughout the decade following tracer addition. In the observations, recovery was also close to 100 % in the days following each tracer addition, but declined to 77 % by the end of the year, largely due to losses during midsummer (Goodale et al., 2015). This N was most likely lost via gaseous pathways (denitrification), as measured leaching losses of ¹⁵N during the first seven months were very low (<1 %). This recovery is very similar to the average recovery rate of 77.9 % for deciduous forests (Templer et al., 2012).

The large difference in total recovery rates between the observations and the simulation originates from the way the ¹⁵N tracer input to the ecosystem is represented in the model. In our simulation the ¹⁵N tracer entered directly into the soluble NO₃ pool of the uppermost soil layer from where it can be taken up by plants or soil microbes, before any gaseous losses occurred. Some gas losses are likely to occur, but because their magnitude cannot be prescribed with confidence, it was not possible to explicitly represent them in the model with this approach. Despite the discrepancy in the total recovery rates, we can use our results to investigate how the ¹⁵N propagates through the different ecosystem compartments.

Because of the way the tracer experiment was represented in the model, the litter layer also did not show the pronounced recovery of the ¹⁵N signal observed in the measurements (Fig. 4). Nevertheless, the model captured several key characteristics of the observations: a large recovery of the signal in the top soil layers and an increasing recovery of the signal with time in the lower soil layers. The simulation supported a long-term retention of N within the ecosystem. This was particularly evident in the scenario results for 2048, as the total recovery rate remained high (88 %). The movement of the ¹⁵N tracer continued into deeper soil layers and the wood became an increasingly important storage pool of the ¹⁵N tracer over time, with only very low losses of the ¹⁵N tracer from the ecosystem. Approximately 11 % of the recovered tracer resided in the plant organs, of which more than 80 % was stored in wood (corresponding to 9 % of total recovery). Due to the high C:N ratio of wood, wood can potentially sequester more carbon per unit of nitrogen than other plant compartments. For NO₃ tracer experiments in temperate forests, Lin et al. (2024) found that almost half of the tracer in plants would reside in stem and branches, that would correspond to wood in our simulations, and more than half in roots and leaves. Our model simulations suggested a larger fraction of the plant compartment's tracer residing in wood than that reported by Lin et al. (2024), implying a greater potential for carbon sequestration. In addition, our simulations cover a longer time span than is currently available from observations.

4.5 Uncertainties in the study

This study has a number of uncertainties. In the measurements the heterogeneity and large mass of the soil pose a considerable challenge (Templer et al., 2012). In the simulation we can capture the overall behaviour based on process understanding. To determine whether the issues identified at Arnot Forest are generic model deficiencies, it would be valuable to extend the analysis to additional experiments.

Since our aim was to assess the N retention at the ecosystem scale and ¹⁵N recovery in the different ecosystem compartments, we consider the model results to be sufficiently robust for this type of analysis. Caldararu et al. (2022) used Latin hypercube sampling to assess the influence of uncertain fractionation rates for different processes (e.g. Robinson, 2001) in the QUINCY model. They found that uncertainties of these fractionation rates did not affect the overall conclusions of their study.

While the QUINCY model covers most fluxes of the terrestrial N cycle, a few processes are missing: it has long been known that some plants are able to also acquire organic nitrogen directly from soil organic matter (Näsholm et al., 1998, 2009), and that symbioses with ectomycorrhizae can also facilitate this uptake (Averill et al., 2014; Chalot and Brun, 1998; Tunlid et al., 2022), but the QUINCY model currently does not include these processes. Their omission might be relevant because half of the trees at the site had ectomycorrhizal associations.

The simulation in this study used the basic soil model (Thum et al., 2019) included in QUINCY. To include abovementioned processes it would be preferable to use as a starting point the more detailed soil model of QUINCY, Jena Soil Model (JSM) (Yu et al., 2020a), which explicitly represents microbial processes and organic matter stabilisation. This new model has been shown to improve the simulation of soil profiles at different sites (Yu et al., 2020a, b) owing to its representation of microbial process and its inclusion of organo-mineral association. Including more explicit process representations instead of relying on first-order kinetics brings the simulation results closer to reality and simulating this tracer study would likely also benefit from using JSM. Using JSM might also help in more realistic representation of immobilization, that has often been identified as a limitation in N cycle modelling studies (Cheng et al., 2019). To fully benefit from the increased process representation of the JSM in future studies, observations of microbial biomass, mineral-associated organic carbon and their ¹⁵N values as well as dissolved organic carbon would be valuable for improving our understanding of the underlying processes.

¹⁵N tracer studies have also shown that the ¹⁵N signal is incorporated into old wood (Tomlinson et al., 2014). This process is not described in the model and would not be captured. In reality, the canopy can also directly take up atmospheric N deposition in some cases (Ferraretto et al., 2022) (but see Da Ros et al., 2023) and thereby influence the carbon uptake of vegetation (Li et al., 2025), even though the soil is the most important pathway for the N entering the plants (Craine et al., 2015). In the Arnot Forest experiment we simulated the ¹⁵N tracer experiment, in which the tracer was sprayed onto the forest floor. The application of the tracer on the forest floor results in lower recovery rates of the tracer than would happen if it was applied on the canopy (Li et al., 2025; Templer et al., 2012). At present we do not have a mechanism in the QUINCY model that would include a mechanism for the uptake of N deposition via leaves.

5 Conclusions

Improving the representation of the N cycle in terrestrial ecosystem models is of high importance, as the N cycle plays a central role in determining future terrestrial C cycle projections. In this study, we simulated a ¹⁵N tracer experiment conducted in a temperate deciduous forest using the QUINCY model, which explicitly represents the ¹⁵N cycle. This enabled us to assess the pathways, turnover rates, and storage of N within the ecosystem. Comparing the simulation with the observations helped to identify priorities for model improvement, including a more detailed evaluation of leaf and fine-root N pools and the vertical transport of N in the soil. Most of the N stored in plant organs was allocated to wood, suggesting a greater C sequestration potential of the trees than if the N had been allocated to other plant organs.

Overall, the model was successful in simulating the main characteristics of the ¹⁵N tracer experiment, taking into consideration some conceptual differences between the model setup and the real-world system. This increases confidence in the model's representation of the N cycle and the ability to simulate the ¹⁵N cycle. QUINCY can therefore be used in studies related to N dynamics and climate change, such as to examine the role of changing N deposition and increasing atmospheric CO₂ over decadal timescales, as studied for example by Bassett et al. (2026). Recent spaceborne observations of δ¹⁵N (Yang et al., 2025) together with QUINCY's ability to simulate the ¹⁵N cycle and leaf chlorophyll, which is also observable from space and simulated by QUINCY (Miinalainen et al., 2025), provide new opportunities to improve our understanding of the global N cycle.

Appendix A

Table A1 Parameters for the calculation of isotopic fractionation from Robinson (2001).