Decadal fates and impacts of nitrogen additions on temperate forest carbon storage: a data–model comparison

To accurately capture the impacts of nitrogen (N) on the land carbon (C) sink in Earth system models, model responses to both N limitation and ecosystem N additions (e.g., from atmospheric N deposition and fertilizer) need to be evaluated. The response of the land C sink to N additions depends on the fate of these additions: that is, how much of the added N is lost from the ecosystem through N loss pathways or recovered and used to increase C storage in plants and soils. Here, we evaluate the C–N dynamics of the latest version of a global land model, the Community Land Model version 5 (CLM5), and how they vary when ecosystems have large N inputs and losses (i.e., an open N cycle) or small N inputs and losses (i.e., a closed N cycle). This comparison allows us to identify potential improvements to CLM5 that would apply to simulated N cycles along the open-to-closed spectrum. We also compare the short- (< 3 years) and longer-term (5–17 years) N fates in CLM5 against observations from 13 long-term ¹⁵N tracer addition experiments at eight temperate forest sites. Simulations using both open and closed N cycles overestimated plant N recovery following N additions. In particular, the model configuration with a closed N cycle simulated that plants acquired more than twice the amount of added N recovered in ¹⁵N tracer studies on short timescales (CLM5: 46±12 %; observations: 18±12 %; mean across sites ±1 standard deviation) and almost twice as much on longer timescales (CLM5: 23±6 %; observations: 13±5 %). Soil N recoveries in simulations with closed N cycles were closer to observations in the short term (CLM5: 40±10 %; observations: 54±22 %) but smaller than observations in the long term (CLM5: 59±15 %; observations: 69±18 %). Simulations with open N cycles estimated similar patterns in plant and soil N recovery, except that soil N recovery was also smaller than observations in the short term. In both open and closed sets of simulations, soil N recoveries in CLM5 occurred from the cycling of N through plants rather than through direct immobilization in the soil, as is often indicated by tracer studies. Although CLM5 greatly overestimated plant N recovery, the simulated increase in C stocks to recovered N was not much larger than estimated by observations, largely because the model's assumed C:N ratio for wood was nearly half that suggested by measurements at the field sites. Overall, results suggest that simulating accurate ecosystem responses to changes in N additions requires increasing soil competition for N relative to plants and examining model assumptions of C:N stoichiometry, which should also improve model estimates of other terrestrial C–N processes and interactions.

1 Introduction

Biogeochemical processes in plants and soils influence Earth's climate by controlling how much carbon dioxide (CO₂) can be removed from the atmosphere and placed into long-term storage in terrestrial ecosystems (Bonan, 2008). Currently, Earth system model ensembles that compare multiple models against each other persistently show a large uncertainty around estimates of CO₂ fluxes exchanged between the land surface and the atmosphere under future scenarios of increasing CO₂ and climate change (Friedlingstein et al., 2006, 2014; Anav et al., 2013). This uncertainty is mainly driven by differences in how models represent biological processes on land and their responses to increasing atmospheric CO₂ concentrations (Lovenduski and Bonan, 2017; Bonan and Doney, 2018). Nutrient limitation is one factor that is likely to constrain CO₂ fertilization of the biosphere (Zaehle and Dalmonech, 2011; Medlyn et al., 2015; Wieder et al., 2015b; Zaehle et al., 2015; Meyerholt and Zaehle, 2018), but not all global land models used in coupled climate–land simulations include explicit representations of the nitrogen (N) cycle. As more global land models add and examine the impacts of coupled carbon (C) and N cycles (Thornton et al., 2007; Sokolov et al., 2008; Zaehle and Friend, 2010; Smith et al., 2014; Goll et al., 2017), additional sources of uncertainty will be added into these models (Wieder et al., 2015a; Lovenduski and Bonan, 2017). Currently, some of the principal uncertainties in simulating terrestrial C–N interactions lie in how models represent plant acquisition of N from soil and the relative competitiveness among plants, decomposers, and denitrifying microbes for soil N (Thomas et al., 2013b, 2015; Medlyn et al., 2015; Zhu et al., 2016b). Thus, evaluating model representations of N cycling is critical for improving our understanding of the magnitude of ecosystem C response to changes in N additions (dC∕dN; Sutton et al., 2008) and how dC∕dN influences the size of the terrestrial C sink over the 21st century.

Human uses of fossil fuels, N-fixing plants, and fertilizers have more than doubled rates of N additions to terrestrial ecosystems compared to preindustrial conditions (Vitousek et al., 1997; Galloway et al., 2003). Increased emissions of reactive N from combustion processes and agricultural sources have led to increases in atmospheric N deposition (Galloway et al., 2003; Vet et al., 2014), which can have multiple effects on forests and other terrestrial ecosystems (e.g., Aber et al., 1998). These effects include shifts in rates of tree growth (Solberg et al., 2009; Thomas et al., 2010) and soil decomposition (Janssens et al., 2010; Frey et al., 2014), as well as increased soil emissions of nitrous oxide (Butterbach-Bahl et al., 2002). The fate of N deposition in plants, soils, or N loss pathways from forests is central to quantifying the effect of N deposition on terrestrial C storage (Emmett et al., 1998a; Nadelhoffer et al., 1999a; Currie et al., 2004; Lu et al., 2010; Templer et al., 2012; Lovett et al., 2013; Wang et al., 2018). Woody plant tissues have higher C:N ratios (e.g., 100–500) than foliage and roots (e.g., 20–40), which allow trees to build more organic C per unit N taken up by plants compared to other plant types. Similarly, woody tissues have C:N ratios that are 1 to 2 orders of magnitude higher than soil organic matter (e.g., 5–25) (Nadelhoffer et al., 1999b; Yang and Luo, 2010; Zechmeister-Boltenstern et al., 2015; Goodale, 2017), allowing trees to store much more additional C if they successfully compete for N deposition than if N is retained in soil; no additional C is stored in forests when N is lost from the system by denitrification or N leaching (Nadelhoffer et al., 1999b).

The fate of N deposition in terrestrial ecosystems has been quantified through field tracer experiments that apply a small amount of highly enriched N with its stable isotope (¹⁵N) to the forest and subsequently measure the recoveries of that ¹⁵N tracer in plant and soil pools. Reviews of these ¹⁵N experiments, which are located predominately in North America and Europe (Tietema et al., 1998; Nadelhoffer et al., 1999b; Curtis et al., 2011; Templer et al., 2012), as well as in warm and humid sites in China (Gurmesa et al., 2016; Wang et al., 2018), indicate that the total amount and partitioning of recovered ¹⁵N varies across sites but that litter and soil pools typically dominate as sinks for N additions during the first few years after tracer application. These litter and soil ¹⁵N sinks often occur directly through microbial or chemical processes within days or weeks after tracer application, without first passing through plants (Berntson and Aber, 2000; Perakis and Hedin, 2001; Providoli et al., 2006; Lewis and Kaye, 2012; Goodale et al., 2015). These ¹⁵N tracer studies are also useful for quantitatively evaluating coupled C–N cycle processes in land models. For example, Thomas et al. (2013b) used mean results from short-term tracer experiments (<3.5 years) to test the responses of two coupled C–N models that treat plant and soil responses to N additions differently: O–CN (Zaehle and Friend, 2010) and the Community Land Model version 4 (CLM4) (Thornton et al., 2007). That analysis showed that CLM4 lost a large fraction of incoming N additions to N gases, while the N retained in the ecosystem was distributed relatively evenly between plants and soils. In contrast, O–CN better estimated total ecosystem retention of N additions and projected that soils dominated the short-term fate of added N.

Comparisons of model simulations with long-term field experiments ultimately provide more relevant constraints on decadal- or centennial-scale forest C uptake and N cycling dynamics than comparisons using short-term field studies that often reflect transient dynamics (Perakis and Hedin, 2001; Jefts et al., 2004; Providoli et al., 2006; Templer et al., 2012). For example, field studies indicate that N initially retained in litter and soil could redistribute to plants and enable additional C uptake over the long-term; alternatively, retained N could accumulate in soil pools or be lost from the ecosystem entirely (Nadelhoffer et al., 2004; Krause et al., 2012; Wessel et al., 2013; Goodale, 2017). However, the long-term fates of N deposition in land models have not yet been evaluated against a synthesis of field measurements. This is in part because, to date, there have only been a handful of individual site-level field studies published that have examined the long-term fates of N additions (e.g., Nadelhoffer et al., 2004; Krause et al., 2012). This study addresses this gap by compiling a summary of ¹⁵N recovery data from long-term ¹⁵N tracer experiments. We then use these data to evaluate the capability of an updated version of the Community Land Model (i.e., CLM5; Lawrence et al., 2019) – the land component of the Community Earth System Model that will be part of the Coupled Model Intercomparison Project Phase 6 (CMIP6) – in its ability to accurately simulate the impacts of ecosystem N additions on annual to decadal timescales. CLM5 includes new, more mechanistic representations of plant N processes, following earlier changes to soil C–N dynamics in CLM4.5 (see Sect. 2.2), which could affect the fate and impact of N in ecosystems. Changes in the fate of N in ecosystems could also be a result of how strongly an ecosystem retains N additions, which is a function of the rates of N added and lost from the ecosystem, as well as internal N cycling dynamics (Aber et al., 1998). To identify whether model-simulated N fates depend on the magnitude of N fluxes, we also compared N fates in CLM5 using two kinds of N cycles – one with high rates of N inputs and losses and one with lower rates of N inputs and losses. This comparison also allowed us to identify potential improvements to the model that are independent of N fluxes. Through this novel data–model comparison project, we provide a synthesis of long-term ecosystem ¹⁵N addition experiments and identify how differences in temporal dynamics of N cycling between field measurements and CLM5 lead to divergences in measured and modeled N fate and ecosystem C responses to N additions.

2 Methods

To assess CLM5's ability to accurately simulate C–N dynamics on both short and longer timescales, we first compiled existing and newly available field data from eight sites that applied a ¹⁵N tracer at least a decade ago. We ran model simulations for each site to examine how N fates and C sink responses to N additions might differ in simulated ecosystems with high N inputs and losses (characteristic of an open N cycle) and simulated ecosystems with low N inputs and losses (characteristics of a closed N cycle). Details on which fluxes the model includes as N inputs and losses are described in Sect. 2.2. Comparing the open and closed representations of the N cycle, described in more detail in Sect. 2.3.1, allows us to examine how sensitive CLM5 is to the “openness” of an ecosystem's N cycle. We also evaluate how the model's C sinks respond to N deposition in both the short term and long term. We define short-term recovery as time points within 3 years after the tracer was applied because the majority of rapid changes in modeled N recovery occur during this time (see Sect. 3) and many ¹⁵N experiments report results within 1–3 years after tracer application (Templer et al., 2012). Longer-term recovery includes time points after 3 (i.e., 5–17) years.

2.1 ¹⁵N tracer field experimental sites

At each of the eight field sites used to evaluate CLM5 (Table 1), a ¹⁵N tracer was added at least 10 years ago, often under both ambient and fertilized conditions. These sites span a range of environmental conditions in North America and Europe and include two plant functional types (PFTs) in CLM5: broadleaf deciduous temperate (BDT) and needleleaf evergreen temperate (NET) trees. Present-day ambient N deposition at these sites ranges from approximately 0.8 to 2.0 g N m⁻² yr⁻¹. Across sites, a ¹⁵N tracer was added as either ammonium, nitrate, or in some cases, in both forms to five ¹⁵N experiments under ambient conditions and to eight experiments under fertilized conditions, with additions ranging from 2.5 to 7.5 g N m⁻² yr⁻¹. Available field measurements of ¹⁵N recovery from these sites are in Table S1 in the Supplement and available at https://doi.org/10.5281/zenodo.2772160.

Table 1 Site information for long-term ¹⁵N experiments in Europe and North America, including corresponding plant functional type (PFT) in CLM5, which is either broadleaf deciduous temperate (BDT) or needleleaf evergreen temperate (NET) tree. The year of stand establishment is the year in which we implemented a harvest in the model to simulate the forest's reported stand age. N deposition is reported as throughfall or the sum of wet and dry N deposition as compiled from available literature. Tracer experiments at plots receiving ambient N deposition and fertilizer treatments are described on separate lines for each site. Please refer to Sect. 2.4 for details on how these field experiments were simulated using CLM5.

^a Magill et al. (2004). ^b Nadelhoffer et al. (2004). ^c Goodale (2017). ^d Elvir et al. (2006). ^e Ivan Fernandez (personal communication, 2018). ^f Nadelhoffer et al. (1999a). ^g Tietema et al. (1998). ^h Gundersen (1998). ⁱ Seftigen et al. (2013). ^j Moldan et al. (2006). ^k Kjønaas and Wright (2007). ^l Krause et al. (2012). ^* We did not simulate the second tracer application that took place in 2012 at the fertilized plots in Bear Brook.

2.2 Model description

We evaluated CLM5 (development version 16_r253) for its ability to estimate the site-level fate of N deposition against the eight experimental sites listed in Table 1. CLM5 is the terrestrial component of the Community Earth System Model (CESM 2.0) and has undergone several changes to its C and N biogeochemistry since CLM4. Briefly, in CLM4.5, the model's original soil biogeochemistry was replaced with a vertically resolved CENTURY-based approach and is described in detail by Koven et al. (2013) and Oleson et al. (2013). In CLM5, three important changes were made to plant C and N dynamics. First, the Leaf Utilization of Nitrogen for Assimilation (LUNA) module allows plants to adjust their photosynthetic capacity (i.e., the maximum rate of carboxylation; V_c,max) based on environmental conditions (Ali et al., 2016). Specifically, V_c,max is influenced by the amount of leaf N allocated for carboxylation, as well as day length and season. Second, plants can alter and optimize their stoichiometry (FlexCN module), which removed the down-regulation of gross primary productivity (GPP) that was used in CLM4 and CLM4.5 (Ghimire et al., 2016). The amount of N that is allocated to individual sub-plant pools is determined based on a fixed set of allometric ratios and the amount of N the plant has for new growth. Additional details on how stoichiometry is optimized can be found in the CLM5 documentation referenced below. Third, in the Fixation and Uptake of Nitrogen (FUN) module, plants pay C costs (which are respired) for acquiring N from symbiotic N fixation, uptake of soil N, and re-translocation (Shi et al., 2015). Additional information about these modifications, as well as other changes to model processes and parameterizations, can be found in the model documentation (Lawrence et al., 2018).

As in prior versions of the model, C and N cycles in CLM5 are coupled at 30 min time steps through plant and soil competition for soil N and internal recycling of plant and soil material through litterfall (Thornton et al., 2007; Koven et al., 2013; Oleson et al., 2013; Thomas et al., 2013a). New N inputs enter an ecosystem through N deposition and free-living and symbiotic N fixation. When N deposition is added to the inorganic soil N pool, it is distributed vertically through the soil column according to an exponential profile; approximately 40 % of N deposition is added to the top 2 cm and approximately 95 % is added to the top 20 cm. When the amount of soil N is smaller than the total N demand, soil N is divided between plants and an implicit representation of microbial immobilization into soil and litter based on each N sink's proportionate demand to the total N demand. Free-living biological N fixation is calculated as a function of annual evapotranspiration and added to the soil mineral N pool. Symbiotic N fixation is passed directly to the plant and depends on plant N demand, the cost of N fixation for the plant, and soil temperature (Lawrence et al., 2018, 2019); details on the model's representation of N fixation are available at https://escomp.github.io/ctsm-docs/doc/build/html/tech_note/FUN/CLM50_Tech_Note_FUN.html (last access: 9 December 2018). Subsequent losses of N occur through production of N gases during nitrification and denitrification. Denitrification occurs in the anoxic portion of the soil and is constrained by decomposition and the availability of nitrate. After gaseous losses, N is lost through water, specifically through surface runoff of dissolved inorganic N over land to stream flow and subsurface leaching through the soil column; the model does not simulate losses of organic N. Rates of decomposition are limited by soil moisture, soil temperature, oxygen availability, and N availability. As litter decomposes into soil organic matter, a portion of C is respired and N is transferred from litter pools through to soil pools. In all our simulations (described below), we turned off transient losses of biomass N from fire and harvest because these disturbances infrequently occur at the sites we simulated.

2.3 Model simulations

For each site and N cycle configuration (see Sect. 2.3.1 below), initial ecosystem C and N stocks for 1850 were generated using a spin-up approach where the model was run using 1850 concentrations of CO₂ (285 ppm) and the model's standard climate forcing dataset from the Global Soil Wetness Project Phase 3 (GSWP3; http://search.diasjp.net/en/dataset/GSWP3_EXP1_Forcing, last access: 6 April 2019) developed by Kim (2017). Specifically, we ran the model for 500 years in accelerated decomposition mode by cycling through the 1901–1920 climate forcing dataset and then for a minimum of 1500 years in regular mode until soil and plant C and N stocks achieved steady state. Subsequently, we ran a historical simulation from 1850 to 2010 (or until 2015 for Arnot Forest, where ¹⁵N recovery was measured after 2010) using transient GSWP3 climate, N deposition, and atmospheric CO₂. We ran CLM5 in single-point mode for each site and modified each site's grid cell area to contain a single PFT and land surface unit (e.g., no lakes). Finally, following Thomas et al. (2013a), we implemented a harvest in the year that established the observed present-day stand age for each site. For Arnot Forest, simulations cycled through the 2006–2010 forcing data for the model years 2011–2015.

2.3.1 Modeling open and closed N cycles in CLM5

Similar to prior versions of the model, CLM5 uses input data for N deposition and models fluxes of other N inputs and losses that are unrealistically high for temperate forests (Thomas et al., 2013a, b). Prior versions of CLM also greatly underestimate measured rates of N losses to leaching and runoff in these ecosystems (e.g., MacDonald et al., 2002; Aber et al., 2003; Nevison et al., 2016; Thomas et al., 2013a). To address the differences between fluxes used in CLM5 and field-based expectations (see Fig. 1), we simulated how each of the field sites in Table 1 would respond to N deposition under (a) CLM5's default N cycle with high N inputs and losses (i.e., open N cycle) and (b) under an alternative “closed” version of the model where N inputs were reduced to better match observations and N losses to denitrification were correspondingly halted to allow the model to build realistic soil N pools and N leaching losses (see Sect. 3.1). In CLM5's configuration with an open N cycle, N deposition rates come from atmospheric modeling simulations described in Lamarque et al. (2010). In the year 1850 (Fig. 1a), the N deposition rate in CLM5 averaged across our sites was 0.4±0.3 g N m⁻² yr⁻¹ (mean ± 1 standard deviation, SD, across modeled sites), exceeding the estimated rate of about 0.1 g N m⁻² yr⁻¹ for preindustrial times (Holland et al., 1999; Fakhraei et al., 2016). The temperate forests we simulated also rarely contain plants with N-fixing symbionts and have rates of free-living N fixation closer to around 0.1 g N m⁻² yr⁻¹ (Vitousek et al., 2013; Tedersoo et al., 2018). However, CLM5 models symbiotic N fixation rates in 1850 of 0.5±0.3 g N m⁻² yr⁻¹ and mean free-living N fixation of 0.3±0.03 g N m⁻² yr⁻¹. In addition, losses of N to denitrification were as high as previously found in CLM4 and CLM4.5 (Thomas et al., 2013a; Houlton et al., 2015), accounting for 99 % of all N losses (Fig. 1b).

Figure 1 Rates of (a) nitrogen (N) inputs and (b) N losses averaged across sites for 1850 (using the last 20 years of the spin-up simulation) and the present day (using the last 20 years of the historical baseline simulation; see Sect. 2.3) for ecosystems with open and closed N cycles. Mean inputs and losses include sources that were >0.1 % of total fluxes. N deposition is based on input data files, while other fluxes are modeled in CLM5 (see Sect. 2.2). For ecosystems using a closed N cycle, N deposition, free-living N fixation, and N deposition were set to observation-based estimates. Denitrification was turned off because of the uncertainty around the portion of N losses due to this loss term (see Sect. 2.3.1). Error bars represent 1 standard deviation simulated across all sites.

In the configuration of CLM5 with a closed N cycle, we changed the model's ecosystem N inputs to be more consistent with historical reconstructions and measurements that suggest 19th century NO_x emissions were smaller (Hoesly et al., 2018) than those estimated by Lamarque et al. (2010). To achieve this, we first lowered preindustrial N deposition to 0.1 g N m⁻² yr⁻¹ for the year 1850 (Fakhraei et al., 2016) and to 0.2 g N m⁻² yr⁻¹ in 1950 to account for the doubling in N deposition and NO_x emissions from widespread use of fertilizers produced by the Haber–Bosch process (Galloway et al., 2003; Reay et al., 2008). Then we set N deposition from 1975 to present to the site-reported N deposition rate (Engardt et al., 2017) in Table 1 and used linear interpolation to calculate N deposition between 1850–1950 and 1950–1975. After 1975, N deposition was held constant to roughly match regional trends (Driscoll et al., 2003; Galloway et al., 2013). These changes allowed us to better simulate the expected N deposition rates at our modeled sites so that the effect of N deposition on N fates between CLM5 and field experiments are more comparable. To adjust N fixation rates, we turned off symbiotic N fixation and set free-living N fixation at all sites to 0.1 g N m⁻² yr⁻¹ to match field expectations (Cleveland et al., 1999; Tedersoo et al., 2018). With these dramatically lower rates for N inputs (Fig. 1a, left panel), simulated C and N stocks were too small, which required us to turn off N losses from denitrification to achieve realistic baseline plant and soil C stocks (Table 3; Sect. 3.1). With denitrification turned off, hydrologic losses of N increased so that the model's present-day subsurface leaching and runoff fluxes each increased to around 0.1 g N m⁻² yr⁻¹ (Fig. 1b, right panel), which is within the range measured in temperate forest stream water (Aber et al., 2003; Gundersen et al., 2006; Groffman et al., 2018). A model that simulates small fluxes of both heterotrophic N fixation inputs and denitrification losses might best match observations of these processes in these temperate forests (Tjepkema, 1979; Roskoski, 1980; Hendrickson, 1990; Barkmann and Schwintzer, 1998; Bernal et al., 2012; Duncan et al., 2013; Morse et al., 2015) but requires substantial model development to achieve (Thomas et al., 2013a; Houlton et al., 2015). Our alternative model included in this study is thus an oversimplification included to examine model C–N responses in ecosystems with a much more closed N cycle than released versions of CLM5, details of which can be found in Lawrence et al. (2018). Our results below highlight how sensitive CLM5 is to the openness of its N cycle, an emergent property that should be the focus of future model development.

2.4 Model N recovery calculations

To calculate N recoveries in ecosystems under ambient N deposition in CLM5, we followed the approach of Thomas et al. (2013b) of adding a small amount of additional N into the N deposition input stream, which enters the model's inorganic soil pool as ${\mathrm{NH}}_{\mathrm{4}}^{+}$. Across the field studies, ¹⁵N tracers were applied differently, with variation in N addition rates, forms (${}^{\mathrm{15}}{\mathrm{NH}}_{\mathrm{4}}^{+}$ or ${}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}^{-}$), and timing (Table 1), although most tracer applications were distributed across the growing season. For our simulations under ambient N deposition, we implemented a consistent approach across all simulations by adding a 0.5 g N m⁻² yr⁻¹ “tracer” in the model during the first year the tracer was applied in the field. In the first year, we applied the N tracer in the model equally across days during April through September to capture the most active portion of the growing season. In CLM5, we ran sensitivity tests for two of our eight model sites (see Table 1 for a full list, including site names), an old growth forest (Alptal) and a younger forest (Harvard NET), which confirmed that the smallest amount of N we could apply while maintaining realistic ecosystem N recovery responses at both sites was 0.5 g N m⁻² yr⁻¹ (Fig. S1 in the Supplement) and is consistent with Thomas et al. (2013b). A sensitivity test for Harvard NET also indicated that the mean N recovery across the last 20 years of the historical simulation was relatively insensitive to which months within the growing season the N tracer was applied (Fig. S2). For simulations with fertilization, we applied the site-reported fertilization rate (Table 1) in the model during all years the fertilizer was applied in the field but only during April through September of each of those years, as we did for simulations under ambient N deposition.

Annual N recovery was then calculated for each N cycle configuration, site, and year by taking the difference in N stocks between a baseline simulation without a tracer or fertilization treatment and its corresponding simulation with tracer or fertilizer added, according to the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathrm{N}\mathrm{recovery}\left(t\right)=\\ & {\displaystyle \frac{\mathrm{N}\mathrm{stock}{\left(t\right)}_{\mathrm{N}\mathrm{addition}}-\mathrm{N}\mathrm{stock}{\left(t\right)}_{\mathrm{baseline}}}{\sum _{t=\mathrm{0}}^t{\mathrm{N}\mathrm{inputs}}_{\mathrm{N}\mathrm{addition}}-\sum _{t=\mathrm{0}}^t{\mathrm{N}\mathrm{inputs}}_{\mathrm{baseline}}}},\end{array}\end{array}$$

where N stock(t)_N addition is the N addition simulation's N stock at year t since the application of the tracer or fertilizer, N stock(t)_baseline is the N stock from the simulation without a tracer or fertilizer added, and N inputs_N addition and N inputs_baseline are the total amounts of N entering the ecosystem in the N addition and baseline simulations at each time, respectively. N inputs are the sum of N deposition, the added tracer or fertilizer, and biological and free-living N fixation. Total vegetation stocks in the model include N in all aboveground and belowground plant pools, including plant stem, foliage, fine roots, and coarse roots. Soil N stocks include N in litter, organic matter, and soil inorganic N pools. We did not include coarse woody debris in the soil stock because coarse woody debris is rarely measured in ¹⁵N tracer experiments. Simulated N recovery in coarse woody debris is possible in CLM5, although this pool accounted for less than 3 % of recovered N. Thus, we applied the field definition of total ecosystem recovery to the model results by defining total ecosystem recovery in the model as the sum of the plant and soil pools. We note that N recoveries of N deposition or fertilizer in ecosystems simulated with open N cycles (reported in Table S1) include a small effect of N fixation rates changing in response to added tracer or fertilizer. Calculated recoveries in simulated ecosystems with closed N cycles (shown in Sect. 3) are a response to the added tracer or fertilizer alone because the control and N addition simulations have the same fixed N fixation rates.

2.5 Calculating ecosystem C response to N additions

To examine the impact of model-estimated N fates on the land C sink, we calculated the modeled change in plant or soil C storage per unit change in N input (g C g⁻¹ N, i.e., dC∕dN), which is frequently done to quantify the impacts of N additions on ecosystem C pools (De Schrijver et al., 2008; Sutton et al., 2008; Thomas et al., 2013a; Frey et al., 2014). Model dC∕dN (g C g⁻¹ N) for each C pool of interest (e.g., total soil or total plant pool) was calculated for each year according to the following equation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{dC}}{\mathrm{dN}}}(t{)}_{\mathrm{model}}={\displaystyle \frac{\mathrm{C}\mathrm{stock}(t{)}_{\mathrm{N}\mathrm{addition}}-\mathrm{C}\mathrm{stock}(t{)}_{\mathrm{baseline}}}{\sum _{t=\mathrm{0}}^t{\mathrm{N}\mathrm{inputs}}_{\mathrm{N}\mathrm{addition}}-\sum _{t=\mathrm{0}}^t{\mathrm{N}\mathrm{inputs}}_{\mathrm{baseline}}}},\end{array}$$

where C stock(t)_N addition is the mean C plant or soil stock at time t in the N addition simulation since the application of the tracer or fertilizer and C stock(t)_baseline is the mean C plant or soil stock at time t in the baseline simulation since the application of the tracer in the N addition simulation. Alternatively, we can estimate both the field and model dC∕dN using the method from Nadelhoffer et al. (1999b):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{dC}}{\mathrm{dN}}}\left(t\right)=\sum _{i=\mathrm{1}}^n{\mathrm{N}\mathrm{recovery}}_i\times {\mathrm{C}:\mathrm{N}}_{\mathrm{pool}\left(i\right)},\end{array}$$

where N recovery is the N recovery of the ith component of a measured plant or soil stock (e.g., foliage, wood, litter), C:N_pool(i) is the C:N ratio of the ith component that makes up a particular stock, and the summation is over i pools that make up the total plant or soil stock. We based our C:N ratios for stocks on site-reported values from the literature (see Table 2). Ecosystem dC∕dN is the sum of plant and soil dC∕dN.

Table 2 Mean C:N values for plant and soil pools reported in site-specific literature and simulated by CLM5 (averaged across sites). Model means (±1 standard deviation across simulated sites) are the average of C:N ratios from the last 20 years of the baseline simulation from simulations with a closed N cycle. Observational means (±1 standard deviation across measured site data) are based off site-reported or field-estimated values for C:N ratios, which are listed for each site, along with their references (which describe the sampling methods for each pool) in Table S3.

^* In CLM5, wood and coarse roots have the same C:N ratios and are split into live and dead pools. In this table, we have listed modeled live biomass as coarse roots and modeled dead biomass as wood.

2.6 Statistical analyses

To determine whether model estimates of ecosystem traits, N recovery, and dC∕dN differ from observations, we used one-sample t tests to identify whether the differences between the model and field values were significantly different from 0. For all statistical analyses, we used the Shapiro–Wilk test to check for normality. When statistical assumptions were not met, we tested for differences using a one-sample Wilcoxon rank-sum test.

3 Results

Below, we first report the effects of N cycle openness (i.e., high or low N inputs and losses) on modeled ecosystem C and N stocks (Sect. 3.1). Next, we compare model output to available observations for (a) the change in N recovery over time at individual sites (Sect. 3.2), (b) the mean change in N recovery in plant and soil pools on short (<3 years) and longer (>3 years; 5 to 17 years) timescales for all tracer experiments in this study (Sect. 3.3.1), (c) the estimated impacts of these N fates on C stocks (i.e., dC∕dN; Sect. 3.3.2), and finally (d) the patterns in N recovery that emerge via PFT and fertilizer treatment (Sect. 3.4).

3.1 N cycle openness: comparison of CLM5 with open and closed N cycles

Ecosystem traits from the end of the historical simulations are reported in Table 3, along with available site measurements. Simulations using a relatively open N cycle had different rates of N inputs and losses in 1850 and throughout the historical simulation than simulations using a relatively closed N cycle (Fig. 1). Observed soil C stocks were typically higher than those modeled under both open and closed N cycles in CLM5 (p<0.01, Table 3). However, simulations with both an open and closed N cycle produced present-day aboveground net primary productivity rates (ANPP), leaf area index, plant C stocks, and plant and soil N stocks that were statistically similar to observations (p>0.05, Table 3).

Table 3 Comparison of mean nitrogen (N) and carbon (C) stocks and annual fluxes in modeled ecosystems with an open and closed N cycle at the end of baseline historical simulations. Sites fall into one of two plant functional types (PFT) in CLM5: broadleaf deciduous temperate (BDT) or needleleaf evergreen temperate (NET) tree. Model values are the mean of the last 20 years of the historical simulation and include total plant (above and belowground plant pools) and soil (excluding coarse woody debris and to the depth of field measurements) stocks. The total depth of the soil column in CLM5 is approximately 7.5 m at our simulated sites. ANPP and C stocks reported in the literature as organic matter or biomass were converted to units of C by assuming that 50 % of the organic matter is C.

NR: not reported in literature. ^a Magill et al. (2004). ^b Zhao et al. (2011). ^c Frey et al. (2014). ^d Nadelhoffer et al. (2004). ^e Magill et al. (1996). ^f Elvir et al. (2006). ^g Nadelhoffer et al. (1999a). ^h Goodale (2017). ⁱ Kjønaas and Stuanes (2008). ^j Emmett et al. (1998b). ^k Per Gundersen (personal communication, 2018). ^l Emmett et al. (1998b). ^m Beier (1998). ⁿ Krause et al. (2012). ^o Patrick Schleppi (personal communication, 2018).

The two model configurations differed by an order of magnitude in their estimated ecosystem N turnover time (i.e., pool size of the ecosystem divided by total N loss fluxes from denitrification, leaching, emissions of N₂O from nitrification, and runoff). Simulations using a closed N cycle had a higher mean ecosystem N turnover time of 6500±5300 years (mean ± 1 SD across simulated sites) compared to 880±370 years in simulations using an open N cycle (Table S2). The wide range in ecosystem N turnover time results from site-specific differences, including variations in factors governing organic N storage in soils (e.g., texture and past disturbance). By comparing the ratio between the mean turnover time of N in the ecosystem to the mean turnover time of N in plants, where plant N turnover time was calculated as the size of the plant N pool divided by the loss of N from the plant pool from litterfall, we quantified how frequently N could cycle through the plant pool. This metric describes the potential for a unit of N to be used by plants to produce C before it is lost from the system. In simulations with a closed N cycle, ecosystem N cycled through plants an average of 210±190 times before it was lost, while it only cycled through plants 26±15 times in simulations with an open N cycle. Consequently, the same unit of N in a closed N cycle has a longer retention time in plants than when the N cycle is open, which could lead the model configuration with a closed N cycle to produce more plant and soil C per unit of N. Given high biases in the historical N deposition, biological N fixation, denitrification, and leaching or runoff rates in the model with an open N cycle (see Sect. 2.3.1), we focus the remainder of the study on results from the model configuration with a closed N cycle. For completeness, results of N recovery for CLM5 using an open N cycle, which is the default configuration of N cycling in CLM5, are in Table S1 and at https://doi.org/10.5281/zenodo.2772160.

3.2 Site-level decadal changes in N recovery

The model–measurement comparisons of N recovery for each site, experimental N treatment, and N cycle configuration are given in Table S1 and shown in Figs. 2 and S3–S6 (data in Table S1 are also available at https://doi.org/10.5281/zenodo.2772160). In this section, we highlight the recovery at Harvard Forest BDT (deciduous) and Harvard Forest NET (evergreen) under ambient N deposition because Harvard Forest provides the longest record for a model–measurement intercomparison with ¹⁵N tracer results. Under a closed N cycle, both Harvard Forest PFTs simulated that plants were the dominant immediate fate of this added N, accounting for 63 % of added N in the BDT stand and 49 % in the NET stand in the first year after N addition. In contrast, field measurements demonstrate that plants acquired <10 % of added tracer and that soils were the dominant sink for ¹⁵N (Fig. 2; Nadelhoffer et al., 2004). After the first year, the model estimated that the N initially taken up by plants moved to soils within 3 years in the deciduous stand and within 5 years in the evergreen stand, after which the recovered N subsequently stayed in the soil pool, accounting for approximately 70 % of the added N at the end of 2 decades. This pattern was typical across the sites we simulated and across ecosystems under open and closed N cycles. However, total ecosystem recovery was generally lower in ecosystems with an open N cycle relative to ecosystems with a closed N cycle (Figs. S3–S6).

Figure 2 Modeled recovery of N additions for Harvard Forest broadleaf deciduous temperate tree (BDT; a) and needleleaf evergreen temperate tree (NET; b). Simulations were done using a closed N cycle and are compared to observations (stacked bars). The modeled soil stock includes the organic soil, inorganic soil, and litter pools (excluding coarse woody debris). Recovery is calculated as the difference in stock size between a control simulation and a simulation with a tracer added as 0.5 g m⁻² between April and September in year 0. Plots of recoveries at all other sites and for the version of CLM5 with an open N cycle at Harvard Forest are shown in Figs. S3–S5.

In our data–model comparison, CLM5 was typically unable to capture inter-site variations in N recovery across PFT and fertilization levels, likely due to both model errors and measurement uncertainty (Figs. 2, S3–S6). For example, the measured total recovery of tracer at Harvard Forest (evergreen forest) appeared to increase with time (Nadelhoffer et al., 2004), which may be a result of changes in the sampling locations of soil cores between sampling events (i.e., soil sampling was done further inside plot boundaries in later years than in early sampling years). Thus, in Sect. 3.3, we discuss changes in the temporal patterns in simulated and observed N recovery averaged across sites, aggregated to the short term (<3 years) and long term (>3 years) to capture the temporal break in slow and fast changes in simulated N recovery. The small number of sites for each forest type (deciduous or evergreen) and fertilizer treatment (ambient or fertilized) also limited statistical comparisons by these factors. Thus, to reach more statistically robust conclusions, we combined both PFTs and fertilizer treatments in Sect. 3.3 and examine qualitative differences by these factors in Sect. 3.4.

3.3 Mean response in CLM5

3.3.1 Change in N recovery in CLM5

Across all sites and treatments, there were 14 field measurements from 13 experiments reporting ¹⁵N tracer recovery within 3 years of the start of N additions (i.e., short term) and 14 field measurements from nine experiments reporting ¹⁵N recovery after 3 years (i.e., long term; Table 1). In the short term, CLM5 with a closed N cycle estimated more than twice the mean plant N recovery of tracer (46±12 %; mean ± 1 SD across simulated sites) than was measured in the field (18 ± 12 %; mean ± 1 SD across field measurements; p<0.001; Fig. 3a). Short-term tracer recovery in soil was modeled to be 40±10 %, compared to the observational mean of 54±22 % (Fig. 3c), and was not statistically different from observations. Tracer recovery in plants decreased over time in both CLM5 and the field studies. In the long term, the closed N cycle led to modeled plant N recoveries (23±6 %) that were closer to observations (13±5 %) than they were on the short term (Fig. 3b) but were still roughly twice the observed values (p<0.001; Fig. 3b). The modeled decrease in mean plant N recovery over time corresponded with an increase in mean soil N recovery (to 59±15 %) that was smaller than observations in the long term (69±18 %; p<0.05; Fig. 3d). The model's initially high plant N recovery and its later increase in soil N recovery indicate that CLM5 estimates soils to become a dominant sink for N in the long term, but this response is a result of an overcompetitive plant pool that transfers recovered N to soils through turnover of plant litter. Similar to the model configuration with a closed N cycle, CLM5 with an open N cycle overestimated short-term (38±7 %; p<0.001) and long-term (20±6 %; p<0.01) plant N recovery and underestimated long-term soil N recovery (p<0.001). However, CLM5 with an open N cycle underestimated (45±15 %) short-term soil N recovery (39±14 %) compared to observations (p<0.05).

Figure 3 Boxplot showing the mean (filled dot), median (horizontal line), 1st quartile, and 3rd quartile of N recovery (%) in the short term (<3 years) and long term (>3 years) from ¹⁵N experiments with available data (Obs) and model simulations using a closed N cycle (Closed CLM5). Field and modeled data are an aggregate of values from both plant functional types (PFTs) and all N fertilization or ambient conditions. Whiskers extend to the minimum and maximum N recoveries that are not outliers, which are represented by open circles. Variation in observations occurs from differences measured across sites; variation in the model data occurs from differences estimated by CLM5 across modeled sites. Different letters indicate groups that are statistically different (p<0.05).

On average, CLM5 with a closed N cycle simulated that most added N remained in the ecosystem over both short and longer timescales (Fig. 3e and f). Within 3 years of simulated tracer or fertilizer addition, the mean whole ecosystem recovery of N under a closed N cycle was 87±14 %, which was higher than the observational mean of 72±23 % (p<0.05). In the longer term, simulations indicate that these forests retained added N with minimal loss (83±17 %), which was similar to the observational mean of 82±16 %. CLM5 with an open N cycle simulated short-term ecosystem recovery (77±14 %) that was similar to observations. However, an open N cycle in CLM5 led to a longer-term ecosystem recovery that was lower than observations (65±16 %, p<0.01).

Figure 4 Mean short-term (<3 years) and longer-term (>3 years) response of C stocks to added N (dC∕dN) estimated from ¹⁵N experiments (Obs) and model simulations using a closed N cycle (Closed CLM5). The dC∕dN for field experiments (leftmost bars in a and b) is calculated by multiplying the site-level ¹⁵N recovery for foliage, wood, bark, fine roots, coarse roots, O horizon, and mineral soil horizon by the C:N ratio for each pool (Eq. 3; see Table S3 for C:N ratios) and then averaging the total plant or total soil dC∕dN across sites. For equivalent comparison to observations, model dC∕dN (center bars in a and b) is calculated indirectly using observational C:N ratios (indirect) and Eq. (3), as well as directly from the model using the model's C:N ratios (direct) according to Eq. (2) (rightmost bars in a and b). For plant dC∕dN, data from Gårdsjön, Aber low, and Aber high experiments were not included because sub-pools from plants were not reported. For soil dC∕dN, data from Gårdsjön and Bear Brook fertilized were not included because sub-pools from soils were not reported. Error bars represent 1 standard deviation. Variation in the observations occurs from differences in N recoveries and C:N ratios measured across plant and soil sub-pools and sites. Variation in the indirect model estimates occurs from differences in N recoveries estimated by CLM5 across modeled sites and from variation in measured C:N ratios. Variation in the direct model estimates occurs from differences in CLM5's stocks and C:N ratios across modeled sites. Statistical differences between observations and model estimates of dC∕dN are indicated with asterisks for plants (above the white bars) and for soils (below the gray bars), where ${}^{**}$ represents p<0.01 and ${}^{***}$ represents p<0.001.

3.3.2 Change in C response to N additions in CLM5

To scale and compare the effect of plant and soil N recoveries on forest C sinks between CLM5 with a closed N cycle and field measurements, we estimated changes in plant, soil, and total C stocks (i.e., sum of plant and soil stocks) in response to N tracer or fertilizer additions: referred to as (dC∕dN)_plant, (dC∕dN)_soil, (dC∕dN)_total, respectively (Fig. 4). For the model, annual dC∕dN values were computed directly (i.e., “direct approach”) as the difference between the total plant or total soil C stocks between the baseline and tracer (or fertilizer) simulation divided by the difference in the amount of cumulatively added N between the two simulations (Eq. 2). To estimate the effect of field-measured ¹⁵N recoveries on forest C pools, we used the scaling exercise presented by Nadelhoffer et al. (1999b) that we described in Eq. (3), where for each experimental time point dC∕dN is estimated for foliage, wood, bark, fine roots, and coarse roots (when available); the O horizon; and the mineral soil, using the measured ¹⁵N recovery in each pool and the published field-measured values of C:N for that site's particular pool (Table 2, Eq. 3). Because differences in model and field estimates of dC∕dN can occur from differences in total N recovery, distribution of recovered N across sub-pools, or C:N ratios of sub-pools, we also used Eq. (3) to compute a second model-based dC∕dN using the same estimates of field-based C:N ratios, except that bark is not modeled in CLM5. This second, indirect approach allows us to remove sub-pool C:N ratios as a confounding factor in estimates of dC∕dN. We used these two methods to calculate model dC∕dN in order to (a) directly show the model's overall C response to N additions and (b) account for substantial differences in modeled and field-based approximations of C:N ratios in plant and soil pools (Table 2). It should be noted that this scaling exercise depends on the accuracy of the C:N ratios measured in the field and operates under the assumption that C:N ratios of plant tissue and soil horizons stay constant over time. Despite these limitations, this budgeting method allows us to roughly compare differences in ecosystem C response to N additions between CLM5 and field measurements.

Averaged across experiments, the short-term direct estimate of (dC∕dN)_plant in CLM5 under a closed N cycle (rightmost white bar in Fig. 4a) was similar to the field-based estimate (leftmost white bar in Fig. 4a), despite the greater than 2-fold difference in plant N recovery between modeled and observed values (Fig. 3a). Within 3 years of N additions, the direct (dC∕dN)_plant in CLM5 was 26±8 g C g⁻¹ N, compared to the field estimate of 19±14 g C g⁻¹ N (Fig. 4a). On the decadal timescale, direct (dC∕dN)_plant in CLM5 became higher than the observational estimate (28±7 and 18±7 g C g⁻¹ N, respectively, p<0.01, Fig. 4b), though not as large as would be expected based on the amount of long-term plant N recovery. Differences between the directly modeled (dC∕dN)_plant and field (dC∕dN)_plant may also be due to differences in C:N ratios of plant sub-pools. In particular, the C:N of wood in CLM5 is substantially lower (∼266–293) than field-based estimates (411±110, Table 2), compensating for the model's overestimate in N recovery. When the same C:N ratios are used for both field- and model-based estimates of (dC∕dN)_plant, the overestimation of N recovery in CLM5 is carried more dramatically into (dC∕dN)_plant: the simulated (dC∕dN)_plant (Fig. 4a and b, center white bars) becomes substantially higher than observations in both the short (52±15 g C g⁻¹ N, p<0.001) and the longer term (45±15 g C g⁻¹ N, p<0.001).

Figure 5 Ratio of plant recovery and soil recovery of added N into (a, c) broadleaf deciduous temperate (BDT) and (b, d) needleleaf evergreen temperate (NET) forests for sites with a closed N cycle under ambient deposition (a, b) and fertilized (c, d) conditions. Circles represent the ratio of plant to soil recovery of ¹⁵N as measured in field experiments. A ratio of 1.0 represents equal recovery of N in plant and soil pools.

In soils, the direct approach to calculating (dC∕dN)_soil (rightmost gray bars in Fig. 4a and b) in CLM5 with a closed N cycle estimated that soil C stocks would decrease within 3 years of N additions ($-\mathrm{2}\pm \mathrm{4}$ C g⁻¹ N) and increase slightly in the long term (5±3 g C g⁻¹ N). This short-term decline in soil C in response to N additions is a result of decreasing litter C stocks. In contrast, observations indicate that soils retain ¹⁵N, which can be associated with C accumulation (12±5 g C g⁻¹ N; Fig. 4a, p<0.001). In the long term, measured soil C stocks (15±5 g C g⁻¹ N; Fig. 4b, p<0.001) increased more than direct estimates from CLM5. When using the same soil C:N ratios to calculate (dC∕dN)_soil in CLM5 as estimated from observations, short- and longer-term soil (dC∕dN)_soil become similar between the model and field estimates (Fig. 4a and b, center bars; 10±3 g C g⁻¹ N and 15±6 g C g⁻¹ N, respectively). Overall, CLM5 directly estimates a short-term (dC∕dN)_total of 24±7 g C g⁻¹ N and a longer-term (dC∕dN)_total of 33±9 g C g⁻¹ N, while the longer-term estimate ranges from 30 to 106 g C g⁻¹ N when using field-estimated C:N ratios. Generally, CLM5 with an open N cycle followed similar patterns as the model configuration with the closed N cycle for short- and long-term plant and soil dC∕dN, except that the indirect calculation of long-term soil dC∕dN was statistically different from observations. Because of existing model limitations in N cycle representation, model-estimated values of dC∕dN are intended to provide a sensitivity test of how the modeling of N fates can affect model estimates of ecosystem C response to N additions relative to what is expected from field measurements.

Table 4 N recovery in CLM5 simulations using an open N cycle (i.e., the default version of CLM5) and a closed N cycle compared to observations for plant and soil stocks in the short (<3 years) and long term (>3 years). Data are separated according to PFT (i.e., broadleaf deciduous temperate, BDT, and needleleaf evergreen temperate, NET, trees) and fertilization treatment, as well as aggregated across all sites (combined across PFTs and fertilization treatment). The number of data points for each N recovery is listed in parentheses.

3.4 Impacts of forest type and fertilization

In the field, forest types might respond to N deposition differently because of differences in their plant and ecosystem traits (Cornelissen, 1996). In CLM5, the evergreen and deciduous PFTs differ, especially in their foliage C:N (see Table 2) and timing of plant N demand, which should alter decomposition and N mineralization. However, statistical comparisons between modeled and measured recoveries of N additions by forest type and fertilizer treatment are difficult to construct because of the small number of sites available for each category. Despite these limitations, we identified a few recurring differences between (a) deciduous and evergreen forests, as well as (b) between ambient N deposition and fertilizer conditions.

Under ambient conditions, simulations using a closed N cycle had plants with notably more mean recovery of added N in two BDT (64 %) forests than in four NET forests (44 %) in the short term (Table 4, Fig. 5). In the long term, the amount of N recovered in modeled plants decreased in both forest types with no difference between the two PFTs (20 %–23 %; Table 4, Fig. 5). Conversely, simulated recovery of added N in soil was higher in NET (52 %, n=4) than in BDT (32 %, n=2) forests in the short term. However, CLM5 estimated similar long-term recoveries of N additions in soil (72 % to 73 %) in both PFTs, similar to long-term patterns of simulated plant N recovery. In contrast, measured recoveries of ¹⁵N did not differ by forest type for plants or soils at either time point, except for short-term soil N recovery (Table 4, Fig. 5a and c). Simulations of CLM5 with an open N cycle followed similar patterns, except there was an underestimation of long-term soil recovery in both PFTs (Table 4).

Fertilization altered the simulated partitioning of N between plants and soils over short and long timescales and generally reduced overall recovery of N in the ecosystem (Table 4, Fig. 5). Simulations with a closed N cycle yielded lower recovery of N in deciduous plants (32 %) than the two unfertilized stands (64 %) in the short term. For the evergreen stands, modeled plant recovery of N did not differ between the four unfertilized and six fertilized stands, although there was considerable variation among the latter. Simulations for both forest types contradict observations showing that fertilization increases short-term plant recovery of tracer N regardless of forest type (Table 4). However, both measurements and simulations using a closed N cycle demonstrated that fertilization led to a decline in the amount of tracer retained in soil in the long term (Table 4). The model configuration with an open N cycle generally followed similar trends as CLM5 with a closed N cycle.

Overall, we find evidence suggesting that (a) under ambient treatment, CLM5 simulates differences in short-term plant N recovery between deciduous and evergreen forests, while measurements show no discernable difference (Table 4), (b) the movement of N from plant to soil pools over the decadal timescale occurs in a distinctly different manner between CLM5 and measurements (Fig. 5), and (c) the model and measurements respond differently to fertilizer in the short term, but both estimate declines in soil N recovery after at least 3 years of fertilizer additions (Table 4).

4 Discussion

This study compares estimates of ecosystem recovery of N deposition between CLM5, a land model with coupled C and N cycles, and long-term ¹⁵N tracer experiments in temperate deciduous and evergreen forests. We examined CLM5, with a focus on simulations with a more closed N cycle, along three important axes of terrestrial C–N modeling: its ability to simulate (a) the decadal patterns of N recovery in plant and soil pools, (b) the plant and soil C responses to the model's estimates of N recovery, and (c) the potential impacts of forest type and increases in N deposition on the partitioning of recovered N in ecosystems. Below, we also discuss the role of total N inputs and losses in ecosystem N recovery. Based on N recovery patterns, we identify some potential causes for the discrepancy between modeled and observed N fates in plants and soils, focusing on plant uptake and soil immobilization processes and recommending changes for modeling plant–soil–microbial competition in future versions of CLM. We then compare our model estimates of the effects of N deposition on forest C sinks with other measurements in the literature and discuss potential mechanisms behind differences in these responses. Lastly, we discuss the ¹⁵N tracer dataset as a tool for evaluating CLM5 and other land models.

4.1 Modeling ecosystem inputs and losses

Our analysis of CLM5 configured with an open N cycle (i.e, the default configuration of N cycling in CLM5) identified that the model continues to have large biases in N losses (Fig. 1), similar to assessments of previous versions of the model (Koven et al., 2013; Thomas et al., 2013b; Houlton et al., 2015; Nevison et al., 2016). Specifically, CLM5 has unrealistically high rates of denitrification and low rates of N leaching and runoff compared to field measurements. We also identified that rates of preindustrial N deposition in the input datasets were higher than expected from reconstructions for the northeastern United States and parts of Europe (Fig. 1; Fakhraei et al., 2016; Holland et al., 1999). Although simulations with an open N cycle shared some similar responses to N additions as simulations with a closed N cycle, having higher N deposition and denitrification fluxes (i.e., the default version of CLM5) typically led to (a) less total ecosystem recovery of N than simulations using a closed N cycle and (b) an underestimation of long-term soil N recovery compared to observations. In adjusting N inputs and losses in CLM5 to better match field expectations, many of the simulated ecosystem stocks and fluxes (i.e., plant N, soil N, plant C, leaf area, and ANPP) remained similar to observations and simulations using an open N cycle (Table 3). Given that the openness of an ecosystem's N cycle changes the ecosystem's recovery of N inputs within a decadal timescale, we suggest that future model development not only test new mechanistic representations of N fixation and losses from denitrification, nitrification, leaching, and runoff but do so in concert with modified N deposition datasets to ensure that both inputs and losses capture field expectations.

4.2 Plant–soil N competition, plant N uptake, and soil N immobilization

Compared to observations, plants in CLM5 are a larger than expected short-term sink for N additions, with soils becoming the dominant sink for N thereafter, as leaf and fine root litter is incorporated into soils. Field experiments demonstrate the opposite pattern, which is that large amounts of N are directly recovered in soils from the start of tracer application, without passing through plants (Emmett et al., 1998a; Gundersen et al., 1998; Tietema et al., 1998; Nadelhoffer et al., 2004; Goodale, 2017). Nadelhoffer et al. (1999b) used earlier ¹⁵N tracer studies to illustrate similar problems in an earlier generation of models 2 decades ago, but issues continue to persist in how these C–N competition processes are represented. The overly strong plant sink for N deposition in CLM5, both when the N cycle is open and closed, likely results in part from how the model handles N competition between plants and soil immobilization, as well as the model's representation of the plant uptake and soil immobilization processes themselves. Our results suggest the need for additional improvements to CLM5's partitioning of N among plants, soils, and N loss pathways, which is similar to results shown in earlier studies with CLM4 (Thomas et al., 2013b), even after substantial changes to the model's soil (Koven et al., 2013) and plant (Lawrence et al., 2018) C–N biogeochemistry have been made since earlier versions.

In CLM5, the amount of N that plants can acquire depends on how much inorganic soil N is available as well as the total demand for N from all modeled ecosystem processes, including soil immobilization, denitrification, and nitrification (Lawrence et al., 2018). When there is not enough N to meet the total demand from both plants and immobilization, inorganic soil N is divided between plants and soils by proportionately scaling their individual demands to the total demand. Plants can then take up their allocated portion of soil inorganic N if they have enough available C to pay for the cost of taking up that N. Given that plants recover too much added N in the short term, regardless of whether the N cycle is open or closed, an option for reducing plant access to N without reducing the availability of inorganic N for immobilization is to increase the costs for plants to acquire N. To date, CLM has also used the long-standing assumption that plants acquire N only from inorganic N pools rather than organic N and that plant demand does not affect N mineralization rates. However, evolving views of plant–soil interactions suggest more complex representations of both processes may be needed, in which plant mycorrhizal associations and priming can enable plants to acquire N from litter and soil organic matter, rather than relying solely on inorganic N (Schimel and Bennett, 2004; Phillips et al., 2013; Tang and Riley, 2014; Terrer et al., 2016; Zhu et al., 2016a; Sulman et al., 2017). Allowing plants to access N from organic as well as inorganic N pools in the model might seem a counterintuitive suggestion, given that plants already show excessive acquisition of newly added inorganic N. But the inclusion of these microbial-driven processes of N acquisition would both better match current understanding of plant–soil–microbial interactions and could allow plants to meet their overall N demand, even if competition for inorganic N by immobilization were to be increased (see below), and allow for added inorganic N to be retained in the soils longer, similar to observations.

Contrary to our model results, field experiments summarized here and elsewhere (e.g., Tietema et al., 1998; Curtis et al., 2011; Templer et al., 2012) demonstrate that soils dominate the fate of added ¹⁵N and that this soil N sink is both rapid and direct without passing through plants. For example, ¹⁵N tracer studies at two sites simulated here, Alptal (Providoli et al., 2006) and Arnot Forest (Goodale et al., 2015), as well as at other temperate forests (Seely et al., 1998; Perakis and Hedin, 2001; Hagedorn et al., 2005; Lewis and Kaye, 2012), show that large quantities of ¹⁵N can be recovered in association with soil organic matter pools within days or weeks of its addition, including in the soil clay or “heavy” fractions, which are generally the most stable components of soil. CLM5 currently immobilizes little N directly into soil, particularly when N input fluxes are high (open vs. closed N cycle and fertilized vs. ambient simulations). However, it should be noted when comparing model estimates and observations to each other that plant access to existing soil N and N additions in the field is not represented identically in CLM5. For example, N additions (e.g., N deposition) in CLM5 are directly added to the dissolved inorganic N pool, which the model immediately distributes throughout the soil column according to an exponential profile. In contrast, the ¹⁵N tracer in the field is typically applied on top of the leaf litter layer. Although the model and field experiments differ in how they apply the N tracer to soil (directly into the inorganic soil N pool vs. to the top of litter, respectively), the large magnitude of the observed soil N sink and the model's poor ability to reproduce it suggests that modeling a stronger soil immobilization sink should be a priority.

Several soil processes and ecosystem traits that are involved in immobilization are not currently represented in CLM5 and could help increase the N demand for soil immobilization. The model's current soil C–N dynamics (Koven et al., 2013) were adapted largely from the CENTURY model, which has implicit microbial processes rather than explicit representation of microbial N uptake and turnover, which are processes that form a dominant pathway for N incorporation into soil organic matter (Bingham and Cotrufo, 2016). Incorporating an explicit representation of microbial biomass and providing microbes with access to inorganic soil N before plants can access it should increase rates of soil immobilization of added inorganic N, particularly because microbial activity and demand for N is greatest at the soil surface (Iversen et al., 2011; Li and Fahey, 2013), where fresh C inputs are greatest, C:N ratios are high, and microbes have the opportunity to rapidly capture N deposition. Previous modeling work has shown that explicitly representing microbes improves soil C stock projections (Wieder et al., 2013) and that more precisely representing plant and microbial biomass and their enzyme affinity for inorganic N better captures the fates of N in grasslands (Zhu et al., 2017).

Field experiments also demonstrate that when ecosystem N additions increase (e.g., with fertilizer), the recovery of ¹⁵N subsequently increases in plants and decreases in soils in the short and long term (Nadelhoffer et al., 1999b; Templer et al., 2012). But under both open and closed N cycles, CLM5 estimated a different short-term plant response depending on forest type, generally estimating a decrease in plant N recovery in deciduous forests and a slight increase in evergreen forests (Table 4). This response is likely due in part to excess plant uptake of N in forests after the tracer is applied in the model. However, the model generally produced a decline in long-term soil N recovery in response to fertilizer treatment (Table 4), except for the mean soil N recovery in evergreen forests with open N cycles. Model development that incorporates plant–soil–microbial dynamics, as described above, would likely yield larger decreases in soil N recovery deposition and fertilizer because soils would recover more N additions under ambient conditions before fertilizer is added. It is important to note that additions of parameters or process-based representation of ecological processes can add uncertainty to model projections. To limit this added uncertainty, new model representations should be designed and evaluated using robust and representative process-based datasets, which are discussed in other modeling papers, including Prentice et al. (2015), Lovenduski and Bonan (2017), Lombardozzi et al. (2018), and Sulman et al. (2018).

4.3 Forest carbon sequestration from N deposition

Despite plants in CLM5 being overly competitive for ecosystem N additions, the model did not dramatically overestimate the response of ecosystem C pools to these additions in simulations with closed N cycles, when compared to observations (Fig. 4a and b). Instead, CLM5 directly estimated a longer term (dC∕dN)_total, ranging from 19 to 45 g C g⁻¹ N in simulations with closed N cycles, close to the lower bound of the 50 to 75 g C g⁻¹ N range of dC∕dN estimated from measurements of forest growth across a N deposition gradient (Sutton et al., 2008) and similar to the range measured in temperate forests (−2 to 48 kg C kg⁻¹ N), summarized by Frey et al. (2014). Similarly, CLM5's estimates of (dC∕dN)_total for simulations with closed N cycles also fall within the range modeled by O–CN (2 to 79 g C g⁻¹ N) and CLM4 (24 to 30 g C g⁻¹ N) for temperate forest ecosystems (Zaehle and Friend, 2010; Thomas et al., 2013a). However, modeled estimates of dC∕dN can be difficult to interpret against field estimates of dC∕dN or those reported by other models because estimates of dC∕dN depend not only on N recovery in plant and soil pools but also on the C:N ratios of these pools.

For example, our directly measured model results for (dC∕dN)_total initially appeared counterintuitive because CLM5 estimated approximately twice as much recovery of N in plants than measured in field experiments. In CLM5, the C:N ratio of deadwood is approximately half the value of what is measured at field sites (Table 2). In addition, plant C:N ratios in the model appear to have an important role in immobilization, as a higher C:N ratio of plant litter in NET forests led to higher rates of immobilization compared to BDT forests (Fig. 2a and b). When the difference in C:N ratios between the model and field measurements were accounted for by using the mean C:N ratios from available data at our sites, the model's estimated (dC∕dN)_total increased to a range of 30 to 106 g C g⁻¹ N, a range that is much larger than observed estimates. While this range overlaps with the observed ranges, the high end of the range exceeds the observations. This discrepancy between observed and modeled C:N ratios would also apply to the default version of CLM5, where ecosystems generally have an open N cycle. Thus, we recommend that options for improving the calculation of C:N ratios of plant pools in CLM5 be explored in order for CLM5 to more accurately and mechanistically model the correct ecosystem C responses to N additions. To accomplish this, additional field measurements would be needed to evaluate changes to model estimates of C:N ratios and to constrain them to reasonable values.

CLM5 simulations with both open and closed N cycles indicated that adding N to temperate forests yielded, on average, a small loss of soil C in the short term (due to declines in the litter pool) and a small increase in the long term (Fig. 4). These dynamics are consistent with modeled relief of N limitation to litter decomposition in the short term (Bonan et al., 2012) and with increased plant net primary production (NPP) in the long term. However, reviews of long-term N addition studies in mature forests show increases in soil C stocks that are associated with reduced rates of decomposition rather than an increase in plant litter production (Janssens et al., 2010; Frey et al., 2014). The increase in soil C stocks might be explained by the changes to the plant–soil–microbial feedbacks described above, in which plant acquisition of N under ambient N availability is mediated by microbial symbionts that drive decomposition, and these processes slow when external N supplies increase. Explicit representation of plant–soil–microbial feedbacks, such as in Sulman et al. (2017), could improve model representation of both soil C responses to N addition as well as the plant and soil N fates discussed above.

4.4 Evaluating N fate in modeled forests

Model evaluation can often be a challenge because of the limitations in the availability and consistency of how field measurements are taken and analyzed. To build a useful dataset that can be leveraged for model evaluation and for potential benchmarking (Luo et al., 2012; Collier et al., 2018), differences between N recovery estimates from different field experiments need to be reconciled. Even when focusing on temperate deciduous and evergreen forests, we found a wide variation in the measurements of N recovery within each forest type. This range in measurements made it difficult to identify how much of the mismatch between measurements and CLM5 was a consequence of model weaknesses alone. For example, calculations of tracer recovery in soils can depend on uncertainties associated with measuring the soil organic N pool, which is usually the largest and most difficult pool to quantify in a terrestrial ecosystem (Nadelhoffer et al., 2004). Accuracy at some sites in part depends on how well-constrained measurements of bulk density are, the number of soil cores taken, and the depth of the soil cores, which ranged across sites from as shallow as 5 cm in Bear Brook Forest to as deep as 50 cm at Arnot Forest. In long-term field studies, changes to field sampling procedures over time add uncertainty to comparisons between observed and modeled temporal trends, as well as uncertainty in evaluating the model's N recovery response to fertilizer and forest type. Increasing the number of long-term ¹⁵N tracer experiments at evergreen and deciduous forests could help constrain estimates of N recovery in plant and soil pools at these two forest types and lead to a more robust dataset for future model evaluations.

In addition, measurement data from other biomes are needed to evaluate the global impact of CLM5's estimates of N recovery in plant and soil pools. Currently, most tracer experiments are in North America and Europe, which represent a subset of ecosystem types and climates that exist globally. One of the few ¹⁵N experiments in tropical forests indicates that soil retention of ¹⁵N is similar to that found in temperate forests (Gurmesa et al., 2016), even though tropical forests are typically limited by phosphorus availability and N (Hedin et al., 2009). Additional tracer experiments in the tropics would allow us to evaluate whether this pattern is anomalous for forests in this biome and whether the response of CLM5's tropical PFTs to N additions is similar to that modeled for temperate PFTs. Finally, increasing the number of sampling events at current and future field sites will expand our capacity to test more nuanced hypotheses about temporal patterns in dC/dN over time, both during the first few years after tracer application when more dramatic changes in ¹⁵N recovery in plant and soil pools occur and for timescales longer than 2 decades.

5 Conclusions

The accuracy of Earth system model projections of land C storage relies on how well land models can simulate the long-term responses of plant and soil C stocks to environmental change, including to shifts in N deposition. To evaluate a land model commonly used in global model intercomparison projects, we simulated temperate and evergreen forests in CLM5 with open and closed N cycles and subsequently compared the modeled fate and effect of N additions on C stocks against measurements from long-term ¹⁵N tracer experiments. Overall, we found that a sizable portion of N additions in CLM5 are taken up first by plants and then moved into long-term soil pools through the recycling of plant litter, which is contrary to field experiments that indicate N deposition is predominately immobilized and retained in soils. Given that CLM5 overestimated plant N recovery in the short and long term, the modeled plant dC∕dN responses were smaller than expected, although the model slightly overestimated long-term plant dC∕dN when a closed N cycle was used. A larger plant dC∕dN did not typically occur because the model's C:N ratios for wood were smaller than those generally found in the field, which compensated for the model's plant pool that was overly competitive for N additions. For similar reasons, CLM5 underestimated soil dC∕dN, even though the modeled mean N recoveries in soil were similar to observations when using a closed N cycle). Overall, our data–model comparison suggests further exploration into (a) more accurate N input data; (b) better representations of N fixation, denitrification, nitrification, leaching, runoff, and C:N ratios in CLM5; and (c) the incorporation of additional plant–soil–microbial processes in land models to increase soil immobilization and reduce initial plant N uptake of N additions.