Carbon–nitrogen interactions in European forests and semi-natural vegetation – Part 1: Fluxes and budgets of carbon, nitrogen and greenhouse gases from ecosystem monitoring and modelling

The impact of atmospheric reactive nitrogen (N_r) deposition on carbon (C) sequestration in soils and biomass of unfertilized, natural, semi-natural and forest ecosystems has been much debated. Many previous results of this dC∕dN response were based on changes in carbon stocks from periodical soil and ecosystem inventories, associated with estimates of N_r deposition obtained from large-scale chemical transport models. This study and a companion paper (Flechard et al., 2020) strive to reduce uncertainties of N effects on C sequestration by linking multi-annual gross and net ecosystem productivity estimates from 40 eddy covariance flux towers across Europe to local measurement-based estimates of dry and wet N_r deposition from a dedicated collocated monitoring network. To identify possible ecological drivers and processes affecting the interplay between C and N_r inputs and losses, these data were also combined with in situ flux measurements of NO, N₂O and CH₄ fluxes; soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching sampling; and results of soil incubation experiments for N and greenhouse gas (GHG) emissions, as well as surveys of available data from online databases and from the literature, together with forest ecosystem (BASFOR) modelling.

Multi-year averages of net ecosystem productivity (NEP) in forests ranged from −70 to 826 g C m⁻² yr⁻¹ at total wet + dry inorganic N_r deposition rates (N_dep) of 0.3 to 4.3 g N m⁻² yr⁻¹ and from −4 to 361 g C m⁻² yr⁻¹ at N_dep rates of 0.1 to 3.1 g N m⁻² yr⁻¹ in short semi-natural vegetation (moorlands, wetlands and unfertilized extensively managed grasslands). The GHG budgets of the forests were strongly dominated by CO₂ exchange, while CH₄ and N₂O exchange comprised a larger proportion of the GHG balance in short semi-natural vegetation. Uncertainties in elemental budgets were much larger for nitrogen than carbon, especially at sites with elevated N_dep where N_r leaching losses were also very large, and compounded by the lack of reliable data on organic nitrogen and N₂ losses by denitrification. Nitrogen losses in the form of NO, N₂O and especially ${\mathrm{NO}}_{\mathrm{3}}^{-}$ were on average 27 % (range 6 %–54 %) of N_dep at sites with N_dep < 1 g N m⁻² yr⁻¹ versus 65 % (range 35 %–85 %) for N_dep > 3 g N m⁻² yr⁻¹. Such large levels of N_r loss likely indicate that different stages of N saturation occurred at a number of sites. The joint analysis of the C and N budgets provided further hints that N saturation could be detected in altered patterns of forest growth. Net ecosystem productivity increased with N_r deposition up to 2–2.5 g N m⁻² yr⁻¹, with large scatter associated with a wide range in carbon sequestration efficiency (CSE, defined as the NEP ∕ GPP ratio). At elevated N_dep levels (> 2.5 g N m⁻² yr⁻¹), where inorganic N_r losses were also increasingly large, NEP levelled off and then decreased. The apparent increase in NEP at low to intermediate N_dep levels was partly the result of geographical cross-correlations between N_dep and climate, indicating that the actual mean dC∕dN response at individual sites was significantly lower than would be suggested by a simple, straightforward regression of NEP vs. N_dep.

1 Introduction

The global terrestrial net sink for atmospheric carbon dioxide (CO₂) is approximately 1.7 Pg C yr⁻¹, i.e. roughly one-fifth of global CO₂-C emissions by fossil fuel combustion and industry (9.4±0.5 Pg C yr⁻¹). This corresponds to the land-based carbon (C) uptake of 3.2±0.8 Pg C yr⁻¹ minus emissions from deforestation and other land-use changes of 1.5±0.7 Pg C yr⁻¹. The ocean sink is of the same order (2.4±0.5 Pg C yr⁻¹), while twice as much CO₂-C (4.7±0.02 Pg C yr⁻¹) is added yearly to the atmosphere (Le Quéré et al., 2018). Data from atmospheric CO₂ inversion methods (e.g. Bousquet et al., 1999; Ciais et al., 2010), from national to global forest C inventory approaches (Goodale et al., 2002; Pan et al., 2011) and from eddy covariance (EC) flux networks (Luyssaert et al., 2007) have suggested that a dominant part of this terrestrial CO₂ sink is currently occurring in forests, and especially in boreal and temperate forests of the Northern Hemisphere (Ciais et al., 2010; Pan et al., 2011). Tropical forest areas are believed to be closer to carbon neutral (Pan et al., 2011), or even a net C source globally (Baccini et al., 2017), due to emissions from deforestation, forest degradation and land-use change offsetting their sink potential. However, others (Stephens et al., 2007) have argued that the tropical land CO₂ sink may be stronger – and the northern hemispheric land CO₂ sink weaker – than was generally believed. At the European scale, Schulze et al. (2010) calculated that the net biome productivity (NBP, the mean long-term carbon sink at a large spatial scale) of temperate and boreal forests was 81 % of the total continental-scale land sink.

The large European and North American CO₂ sinks have been attributed to a combination of factors including afforestation of abandoned land and formerly cut forests, reduced forest harvest, CO₂ fertilization, changes in management and age structure legacy effects in Europe (Vilén et al., 2016), and atmospheric reactive nitrogen (N_r) deposition (Reay et al., 2008; Ciais et al., 2013, and references therein; De Vries et al., 2017). However, some studies (Nadelhoffer et al., 1999; Gundale et al., 2014; Fernández-Martínez et al., 2017) have questioned the widespread theory that elevated N_r deposition boosts forest C sequestration, and the magnitude of the N fertilization effect on forest C sequestration has been a matter of much debate (Magnani et al., 2007, 2008; Högberg, 2007; De Schrijver et al., 2008; de Vries et al., 2008; Sutton et al., 2008; Dezi et al., 2010; Binkley and Högberg, 2016). A better understanding of the impact of nitrogen deposition on natural and semi-natural ecosystems, in particular over forests, and the impact on the carbon and nitrogen cycles as an indirect effect resulting from anthropogenic activities (Canadell et al., 2007) remains key to improving the forecast of regional (de Vries et al., 2017) and global (Du and de Vries, 2018) models.

The relevance of N_r deposition for the global C sequestration potential, or more explicitly the dC∕dN response (change in C storage with change in N_r deposition), has been estimated typically through meta-analyses of N_r addition experiments (e.g. Schulte-Uebbing and de Vries, 2018), or by combining forest growth inventories, together with estimates of Nr deposition obtained from large-scale forest monitoring plots (Solberg et al., 2009; Laubhann et al., 2009; De Vries et al., 2008). Both methods have many sources of uncertainty. One key difficulty in the latter approach lies in estimating total (wet + dry) N_r deposition (N_dep), especially dry deposition, which is highly variable spatially, very challenging to measure and consequently hard to parameterize in regional-scale chemical transport models (CTMs) (Flechard et al., 2011; Simpson et al., 2014; Schwede et al., 2018). The annual or long-term dry deposition component of N_dep to forests, in all the diversity of N-containing forms (gaseous vs. aerosol, reduced vs. oxidized, inorganic vs. organic, e.g. Zhang et al., 2009), has been actually measured (by micrometeorological methods) in very few forests worldwide (Neirynck et al., 2007; Erisman et al., 1996). Due to the large diversity of atmospheric compounds that contribute to total N_r and the complexity of the measurement techniques required for each compound (Flechard et al., 2011), it is even debatable that complete measurements of all N_r deposition terms have ever been achieved anywhere. Thus virtually all studies of the forest dC∕dN response so far have relied on modelled atmospheric N_r deposition estimates, at least for the dry and occult deposition fractions, and further the N_r deposition data being used were systematically provided by the outputs of large-scale regional (e.g. Sutton et al., 2008; Fernández-Martínez et al., 2017) or even global (Fleischer et al., 2013) models, with resolutions of typically 10 km × 10 km or 1^∘ × 1^∘, respectively. Grid averaging in such large-scale models introduces a large uncertainty in local (ecosystem-scale) N_r dry deposition rates (Schwede et al., 2018), particularly when the forest sites are located near agricultural or industrial N_r sources (Loubet et al., 2009; Fowler et al., 1998).

Additionally, nitrogen losses may significantly offset atmospheric N_r inputs at eutrophicated and acidified sites, with the consequence that dC∕dN may correlate better with net, rather than gross, atmospheric N_r inputs. Depending especially on the extent of ecosystem N saturation (De Schrijver et al., 2008), substantial N losses may occur in the form of nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) leaching (Dise et al., 2009), nitric oxide (NO) and nitrous oxide (N₂O) emissions (Pilegaard et al., 2006); ammonia (NH₃) bidirectional exchange (Hansen et al., 2013); and emissions of di-nitrogen (N₂) from total denitrification (Butterbach-Bahl et al., 2002) (Fig. 1). The implication is that the carbon response to N_dep would be non-linear, with larger dC∕dN at low N_dep rates and a lowering of dC∕dN as N_dep increases, as suggested in the review by Butterbach-Bahl and Gundersen (2011) and further elaborated in De Vries et al. (2014). The latter authors show in their review that above a certain N deposition level, the dC∕dN response declines due to adverse effects of excess N_r deposition and high soil ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) concentration and nitrification (e.g. acidification, nutrient base cation losses, aluminium mobility), which are known to reduce soil fertility and affect ecosystem health and functioning (Aber, 1992).

Figure 1 Flux terms and boundaries of the carbon (a) and nitrogen (b) budgets discussed in this paper. Net ecosystem productivity NEP = GPP − R_eco (≈ NPP − R_het) based on multi-annual eddy covariance CO₂ flux data. The net ecosystem carbon balance (NECB) includes in addition other C loss fluxes such as DIC/DOC, CH₄ and VOC, as well as harvest, thinning or other disturbances (e.g. fire). Inorganic reactive nitrogen (N_r) budget = N_dep − DIN_leach − NO − N₂O. The total N budget includes in addition organic nitrogen deposition (WSON) and leaching (DON), as well as N₂ inputs and losses from biological fixation and denitrification, respectively. CLBS, CSOM, CR, CLITT: carbon stocks in leaves, branches and stems; in soil organic matter; in roots; and in litter layers, respectively. Terms highlighted in red indicate that direct or measurement-based estimates were not available for some or all sites in our datasets (see also Table 2 for a list of acronyms, Table 3 for a summary of methods and Table S6 for data availability).

Carbon losses through dissolved organic carbon (DOC) and biogenic dissolved inorganic carbon (DIC) leaching can also be significant, especially for wetlands (Dinsmore et al., 2010) and also grassland and cropland ecosystems (Kindler et al., 2011; Gielen et al., 2011). This is relevant for the net ecosystem carbon balance (NECB) or the net biome productivity (NBP) estimates obtained on the basis of EC flux systems and needs to be accounted for as a part of the net ecosystem productivity (NEP) that is not actually stored in the system (Chapin et al., 2006; Schulze et al., 2010) (Fig. 1). Dissolved and/or emitted methane (CH₄) may further represent a significant loss from organic soils (Hendriks et al., 2007), while CH₄ oxidation, which is often observed in well-aerated soils and can be suppressed by N_r addition, especially ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Steudler et al., 1989), may affect the net greenhouse gas (GHG) budget. Nitrogen-deposition-induced N₂O emissions from the forest floor (Pilegaard et al., 2006; Liu and Greaver, 2009), or from denitrification triggered by deposited ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in peatland (Francez et al., 2011), can also offset the gain in the ecosystem GHG balance resulting from a hypothetical nitrogen fertilization effect.

Nitrogen deposition or addition is known to affect soil microbial C cycling in many different ways, for example high-level N enrichment generally leading to reduced microbial biomass and suppressed soil CO₂ respiration (Treseder, 2008); a reduction of basal respiration without significant decline in total microbial biomass, following N addition to incubated peat cores (Francez et al., 2011); and added ${\mathrm{NO}}_{\mathrm{3}}^{-}$ altering directly the oxidative enzyme production by microbial communities and hence controlling extracellular enzyme activity (Waldrop and Zak, 2006). Nitrate addition can lead to a reduction in CH₄ emissions from wetlands and peatlands (Francez et al., 2011), since in anaerobic conditions and in the presence of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as an electron acceptor, denitrifying bacteria can oxidize organic C substrates (e.g. acetate) and thus outcompete methanogenic communities (Boone, 1991). However, if chronic N enrichment of peatland ecosystems leads to floristic changes, especially an increase in vascular plants at the expense of bryophytes, the net effect may be an increase in CH₄ emissions (Nykänen et al., 2002), as the aerenchyma of tracheophytes provides a direct diffusion path to the atmosphere for soil-produced CH₄, bypassing oxidation in the peat by methanotrophs. Excess-nitrogen-induced vegetation composition changes in Sphagnum moss peatland are believed to reduce C sequestration potentials, and the effect is likely to be exacerbated by climate change (Limpens et al., 2011).

This complex web of interactions between the C and N cycles and losses shows the need for integrated approaches for studying the impacts of N_r deposition on C sequestration and net GHG budgets. Ideally, all C and N gain and loss pathways (including infrequently or rarely measured fluxes such as N_r dry deposition, organic C and N leaching fluxes, and GHG fluxes; see Fig. 1) should be quantified at long-term experimental sites to improve and calibrate process-based models. Closing the C and N budgets experimentally at each site of large (e.g. FLUXNET) monitoring networks is unlikely to occur in the near future, but realistic and cost-effective measurement approaches can be used to progressively reduce the uncertainties for the large terms of the budgets. Such approaches were tested and implemented in this study, as part of a large-scale effort, within the NitroEurope Integrated Project (NEU, 2013; Sutton and Reis, 2011), to quantify N_r deposition and N losses from ecosystems, in parallel and coordinated with the CarboEurope Integrated Project (CEIP, 2011) to estimate the net C and GHG balance, for forest and semi-natural ecosystems in Europe.

A main objective of this paper is to build tentative C, N and GHG budgets, as well as analyse C–N interactions empirically, for a wide range of European monitoring sites, by using measurements or observation-based data as far as possible, complemented by modelling. Important methodological goals are to critically examine uncertainties in measurement methods and elemental budgets, to identify knowledge and data gaps, and to assess the current state of process understanding as encoded in models. To this end, we compiled the C, N and GHG flux data from NEU, CEIP and other complementary datasets, using a combination of in situ measurements, empirical relationships, ecosystem modelling, literature and database surveys, at the scale of the CEIP and NEU flux monitoring networks. This study presents the methodologies and discusses the different terms of the budgets, including atmospheric deposition from gas, aerosol and precipitation N_r concentration monitoring, soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching measurements and modelling, GHG and N_r emission estimates from chamber measurements and laboratory-based soil bioassays, EC-tower-based C budgets, and historical published data. Forest ecosystem modelling (BASFOR) is used to simulate C, N and GHG fluxes, with the double objective to compare with actual measurements and to fill some gaps in the datasets. Wherever possible, alternative measurements, datasets or modelled data are shown alongside the primary data in order to provide an estimate of the uncertainty in the different terms. In the companion paper (Flechard et al., 2020), the response of C sequestration to N_dep is quantified using the same datasets.

2 Materials and methods

2.1 Monitoring sites

The study comprised 40 terrestrial ecosystem-scale, carbon and nitrogen flux monitoring sites, including 31 forests (F) and nine natural or semi-natural (SN) short vegetation ecosystems, primarily moorlands, wetlands and extensively managed, unfertilized grasslands (Table 1). The sites spanned a European geographical and climatic gradient from the Mediterranean to the Arctic and from the Atlantic to western Russia (Fig. S1 in the Supplement), an elevation range of −2 m to 1765 m a.m.s.l., a mean annual temperature (MAT) range of −1.0 to 17.6 ^∘C, and a mean annual precipitation (MAP) range of 500 to 1365 mm. Selected references are provided for each site in Table S1 in the Supplement. A list of the main acronyms and abbreviations used in the paper is provided in Table 2.

Table 1 Overview of ecosystem and climatic characteristics and inter-annual mean ecosystem–atmosphere exchange fluxes for forest and semi-natural short vegetation sites.

^a PFT (plant functional types): DB – deciduous broadleaf forest; EN – evergreen needleleaf coniferous forest; EB – evergreen broadleaf Mediterranean forest; MF – mixed deciduous–coniferous forest; SN – short semi-natural, including moorland, peatland, shrubland and unimproved/upland grassland. ^b Maximum canopy height. ^c Maximum leaf area index, defined as one-sided or half of the total. ^d Above mean sea level. ^e Mean annual temperature. ^f Mean annual precipitation. ^g Nitrogen deposition. ^h Gross primary productivity. ⁱ Ecosystem respiration. ^j Net ecosystem productivity. n/a: not applicable.

The forest sites of the study ranged from very young (< 10 years old) to mature (> 150 years old) and can be broadly classified into four plant functional types (PFTs) or five dominant tree categories (Table 1): deciduous broadleaf (DB), evergreen needleleaf (EN, comprising mostly spruce and pine species), mixed deciduous–coniferous (MF), and Mediterranean evergreen broadleaf (EB). Forest species composition, stand characteristics, C and N contents of different ecosystem compartments (leaves, wood, soil), soil physical properties and micro-climatological characteristics are described in Tables S2–S5. Semi-natural short vegetation ecosystems included unimproved (mountainous and semi-arid) grasslands, wetlands and peatlands; they are included in the study as unfertilized, C-rich soil systems, providing a contrast with forests where storage also occurs above ground (thus with different C∕N ratios). Among the 40 EC-CO₂ flux measurement stations, most sites (36) were part of the CEIP CO₂ flux network. A further three CO₂ flux sites were operated as part of the NEU network (EN2, EN16, and SN3), and one site (DB4) was included from the French F-ORE-T observation network (F-ORE-T, 2012). Table S6 provides an overview of the available C, N and GHG flux measurements, detailed hereafter.

Table 2 Main acronyms and abbreviations used in the study.

2.2 Nitrogen fluxes

Input and output fluxes of the ecosystem nitrogen and carbon budgets are represented schematically in Fig. 1. The following sections describe the methods used to quantify the different terms, summarized in Table 3.

2.2.1 Atmospheric deposition

To obtain realistic estimates of total (dry + wet) N_r deposition at the 40 sites of the network, it was necessary to measure ambient air concentrations of the main N-containing chemical species at each location, due to the large spatial heterogeneity in gas-phase concentrations, especially for NH₃. The requirement for local measurements of wet deposition was relaxed because this is much less spatially variable. For both dry and wet components, measurements had to be complemented by models, either to calculate fluxes based on local concentration data at each site or to obtain local estimates from a large-scale CTM when data were missing.

Table 3 Summary of the main methods used to quantify carbon, nitrogen, and greenhouse gas fluxes and budgets for the 31 forests and nine short semi-natural vegetation sites included in this study. Horizontal bars (green: forests; blue: short semi-natural vegetation) indicate the percentage of study sites with available data (filled bars) or without available data (open bars). See also Tables S6–S7 for details at individual sites.

¹ Aubinet et al. (2000). ² Dämmgen (2006). ³ Dinsmore et al. (2010). ⁴ Dise et al. (2009). ⁵ Flechard et al. (2011). ⁶ Gielen et al. (2011). ⁷ Hendriks et al. (2007). ⁸ Ilvesniemi et al. (2009). ⁹ Kindler et al. (2011). ¹⁰ Kowalska et al. (2013). ¹¹ Legout et al. (2016). ¹² Luo et al. (2012). ¹³ Pilegaard et al. (2006). ¹⁴ REddyProc (2019). ¹⁵ Schaufler et al. (2010). ¹⁶ Simpson et al. (2012). ¹⁷ Tang et al. (2009). ¹⁸ van Oijen et al. (2005). ¹⁹ See Table S7.

Atmospheric inorganic N_r concentrations, available from the NEU (2013) database, were measured monthly for 2–4 years in the gas phase (NH₃, HNO₃, HONO) and in the aerosol phase (${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$), using DEnuder for Long-Term Atmospheric sampling (DELTA) systems (Sutton et al., 2001; Tang et al., 2009). Concentrations of nitrogen dioxide (NO₂), not covered by DELTA sampling, were measured by chemiluminescence at a few sites only and were otherwise taken from gridded concentration outputs of the European-scale EMEP CTM (details given below). The N_r data initially reported in Flechard et al. (2011) covered the first 2 years of the NEU project (2007–2008); here, the data from the entire 4-year NEU monitoring period (2007–2010) were used and averaged to provide a more robust long-term 4-year estimate of N_r dry deposition. The inferential modelling method was used to calculate dry deposition for N-containing gas and aerosol species, whereby measured ambient N_r concentrations were multiplied by a vegetation-, meteorology- and chemical-species-dependent deposition velocity (V_d) (Flechard et al., 2011, 2013; Bertolini et al., 2016; Thimonier et al., 2018). In the case of NH₃, a canopy compensation point scheme was applied in some models, allowing bidirectional exchange between the surface and the atmosphere. Considering notoriously large uncertainties in deposition velocities and large discrepancies between the surface exchange schemes currently used in different CTMs, we tried here to minimize such uncertainties by using the ensemble average dry deposition predicted by four different models, as in Flechard et al. (2011).

The dry deposition of atmospheric organic N_r (ON) species not accounted for by the EMEP model (e.g. amines, urea), and not included in DELTA measurements, can contribute a fraction of total N_r deposition. However, Kanakidou et al. (2016) suggest that particulate ON largely dominates the atmospheric ON load, and for particles the main atmospheric removal mechanism is through precipitation. Thus, dry deposition of ON is expected to be much smaller than wet deposition of water-soluble organic compounds (see below).

For wet deposition, several sources of data were used, and the final wet deposition estimate was derived from the arithmetic mean of the different sources, where available. First, within the NEU project, a survey was made of the available national and/or transnational (e.g. EMEP, 2013; ICP, 2019) wet deposition monitoring network concentration data for inorganic N (${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$) in the different European countries hosting one or several CEIP/NEU flux sites. These data were checked for consistency and outliers, harmonized, and then spatially interpolated by kriging to provide measurement-based estimates of solute concentrations in rainfall for each of the 40 sites of this study. Wet deposition was then calculated as the product of interpolated concentration times measured precipitation at each site.

Next, 13 sites (DB1, DB3, DB4, EN4, EN9, EN13, EN14, EB2, EB3, MF1, MF2, SN3, SN8) were identified as lacking local or nearby wet deposition measurements. These sites were equipped for three years (2008–2010) with bulk (open funnel) precipitation samplers (Model B, Rotenkamp, Germany; Dämmgen, 2006), mounted above the canopy or inside a clearing for some of the forest sites, with monthly sample change and analysis. The precipitation samples were stabilized by addition of thymol at the beginning of each exposure period and were analysed subsequently for inorganic N_r (${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$) as well as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, Cl⁻, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, base cations (Mg²⁺, Ca²⁺, K⁺, Na⁺) and pH. A few other sites (EN2, EN8, EN10, EN16, DB2, SN9) were already equipped with wet-only or bulk precipitation collectors. No correction was applied to the bulk deposition estimates to account for a possible contribution by dry deposition within the sampler glass funnel (e.g. Dämmgen et al., 2005), since there did not appear to be any systematic overestimation compared with wet deposition estimates from the monitoring networks or EMEP data (see Results and Fig. S2), even if a more significant bias may be expected in dry (Mediterranean) regions.

In addition to inorganic nitrogen, the wet deposition of water-soluble organic N_r (WSON) compounds was also investigated in precipitation samples at 16 sites (Cape et al., 2012). However, since WSON data were not available for all sites and the measurements were subject to considerable uncertainties (Cape et al., 2012), and also because the contribution of WSON to total N_r deposition was on average less than 5 %, WSON was not included in the final estimates of total N_r deposition.

The last data source was the ca. 50 km × 50 km gridded modelled wet inorganic N_r deposition (also NO₂ concentrations, discussed above), simulated by the European-scale EMEP CTM (Simpson et al., 2006a, b, 2012, 2014) for the years 2007–2010, available from EMEP (2013). The data were downloaded in 2013, and it should be noted that in this data series different model versions were used for the different years. This leads to some uncertainty, especially in the dry deposition estimates, but it is hard to say which model version is the most realistic. Evaluation of the model against measurements over this period has shown quite consistent results for the wet-deposited components and NO₂ concentrations, but the dry deposition rates cannot be evaluated versus actual measurements at the European scale. We chose therefore to make use of all versions and years, giving a small ensemble of simulations.

2.2.2 Soil gaseous and leaching losses

Nitrogen losses to the atmosphere (gaseous emissions) and to groundwater (N leaching) are especially hard to quantify and thus typically cause large uncertainties in ecosystem N budgets. These N_r losses were estimated by direct flux measurements or by indirect empirical methods. Soil NO and N₂O emissions were measured in the field using closed static and dynamic chamber methods, as part of NEU (e.g. EN2, EN10, EN16, DB2, SN3, SN8, SN9) and/or collected from the literature (e.g. EN2, EN10, EN14, EN16, DB2, Pilegaard et al., 2006; long-term data at EN2 in Luo et al., 2012). Such data were available for N₂O at seven forest sites and four semi-natural sites, as well as at five forest sites for NO (Table S6). Manual static chamber N₂O measurements were made manually at a typically fortnightly (growing season) or monthly (winter half-year) frequency at many sites. Automatic chamber systems, allowing continuous N₂O measurements at a frequency of four times per day, were deployed at EN2, EN10, DB2 and SN3. Fluxes of NO were only measured by automatic dynamic (open) chambers. Measured fluxes were scaled up to yearly values by linear interpolation or using the arithmetic mean of all flux measurements. There may be considerable uncertainty in the annual flux if gap-filling is based on linear interpolation between discrete values, when flux measurements are made manually and are therefore discontinuous and infrequent (Parkin, 2008). This is due to the episodic nature and log-normal distribution of NO and N₂O emissions, observed particularly in fertilized croplands and grasslands. However, this episodicity is less pronounced in semi-natural ecosystems, or at least the magnitude of the episodic fluxes is generally much smaller than in fertilized agro-systems (Barton et al., 2015). The uncertainty in annual emissions estimated in our study from manual chamber measurements is related to the observation frequency (fortnightly or monthly) and is likely larger than in the case of automatic (continuous) chamber measurements.

Direct in situ N_r and non-CO₂ GHG gas flux measurements were unavailable at many sites. These soil N₂O and NO (and also CH₄) fluxes were therefore also estimated, as part of NEU, from empirical temperature and moisture responses of soils. These responses were established in a series of factorial soil incubation experiments in controlled conditions with four levels of temperature (5–20 ^∘C) and water-filled pore space (20–80 WFPS %), following the protocol described in Schaufler et al. (2010). Twenty-four undisturbed soil cores (top 5 cm of the mineral soil, Ah horizon) were taken from each of 27 forests and eight semi-natural sites in spring after soils had warmed up above 8 ^∘C for 1 week in order to guarantee phenological comparability of the different climatic zones. Sampling was conducted in 2008, 2009 and 2010, and cores were sent to a common laboratory at the Federal Research and Training Centre for Forests (BFW, Vienna, Austria) for the controlled environment bioassays, which were carried out straight away. The 5 cm topsoil layer was selected as it represents the highest microbial activity and correspondingly high GHG production/consumption rates, although processes in deeper soil layers should not be neglected (Schaufler et al., 2010). Site-specific, empirical bivariate (temperature, WFPS) relationships describing soil fluxes for CO₂, N₂O, NO and CH₄ were derived from the incubation results and then applied to multi-annual time series of soil temperature and moisture measured at the sites, mimicking field conditions and providing scaled-up estimates of potential annual trace gas emissions.

Leaching of dissolved inorganic nitrogen (DIN = ${\mathrm{NH}}_{\mathrm{4}}^{+}+{\mathrm{NO}}_{\mathrm{3}}^{-}$) was measured using lysimeter setups, or estimated from a combination of suction cup measurements (typically ∼1 m soil depth) and a hydrological drainage model, at a few sites during the NEU monitoring period (EN2, EN4, EN10, EN15, EN16, DB1, DB2) and as part of parallel projects (EN8, DB4). One-dimensional (1-D) drainage models were based on the soil water balance equation using evapotranspiration, observed precipitation and changes in soil water content (Kindler et al., 2011; Gielen et al., 2011). For the forest sites where no leaching measurements were available, the empirical algorithm by Dise et al. (2009) was applied to predict DIN leaching based on key variables (throughfall inorganic N_r deposition DIN_TF, organic horizon C∕N ratios, MAT). The algorithm, developed from the extensive Indicators of Forest Ecosystem Functioning (IFEF) database (> 300 European forest sites), simulates the non-linearity of DIN leaching with respect to DIN_TF and soil C∕N ratio, with critical thresholds for the onset of leaching of DIN_TF=0.8 g N m⁻² yr⁻¹ and C∕N = 23, respectively. Since the algorithm requires DIN_TF as input, as opposed to total (above canopy) N_dep, in the present study we applied a reduction factor of 0.85 from N_dep to DIN_TF (i.e. a canopy retention of 15 % of atmospheric N), which was calculated as the average of all available individual DIN_TF ∕ N_dep ratios in the IFEF database. A comparison with values of DIN_TF ∕ N_dep ratios actually measured at the EN2, EN8, EN10, EN16 and DB2 sites (0.71, 0.80, 0.29, 0.85 and 1.11, respectively; mean ± SD 0.75±0.30) shows that the applied ratio of 0.85 is plausible but also that much variability in canopy retention/leaching may be expected between sites.

2.3 Carbon fluxes

2.3.1 Ecosystem–atmosphere CO₂ exchange

Half-hourly rates of net ecosystem–atmosphere CO₂ exchange (NEE) were measured over several years (on average 5 years; see Table S6) by the eddy covariance (EC) technique at all sites. The long-term net ecosystem productivity is defined following Chapin et al. (2006) as the difference between gross primary production (GPP) and ecosystem respiration (R_eco) and is thus calculated as the straightforward annual sum of NEE fluxes (with opposite sign). The net ecosystem carbon balance may differ from the NEP if C fluxes other than assimilation and respiration, such as DIC/DOC leaching, CH₄ and other volatile organic compound (VOC) emissions, as well as lateral fluxes (harvest, thinning) and other disturbances (fire), are significant over the long term (Chapin et al., 2006). For convenience in this paper, we use the following sign convention for CO₂ fluxes: GPP and R_eco are both positive, while NEP is positive for a net sink (a C gain from an ecosystem perspective) and negative for a net source. Previous studies have normalized C flux data through the carbon use efficiency (CUE), commonly defined from a plant's perspective as the ratio of net to gross primary productivity (NPP ∕ GPP), or the biomass production efficiency (BPE = BP ∕ GPP; Vicca et al., 2012), which is a CUE proxy. By analogy, we define here an ecosystem-scale, medium-term indicator of carbon sequestration efficiency (CSE) as the NEP ∕ GPP ratio, calculated from measurable fluxes over the CEIP/NEU project observation periods.

The EC technique is based on fast-response (sampling rates typically 10–20 Hz) open-path or closed-path infrared gas analyser (IRGA) measurements of turbulent fluctuations in CO₂ concentration (c) in the surface layer above the ecosystem, coupled with ultrasonic anemometer measurements of the three components of wind (u, v, w) and temperature. The NEE flux is calculated as the average product of c and w fluctuations, i.e. the covariance (Swinbank, 1951; Lee et al., 2004).

The EC-CO₂ flux measurements reported here followed the protocols established during the CEIP project, largely based on the EUROFLUX methodology (Aubinet et al., 2000). Briefly, post-processing of the raw high-frequency EC data included typically de-spiking to remove outliers, 2-D rotation of the coordinate system, time lag optimization by maximization of the covariance between CO₂ concentration and the vertical component of wind speed (w), and block-averaging over the flux-averaging interval of 30 min. Corrections were applied for various methodological artefacts, including notably (i) flux losses at the different frequencies of flux-carrying eddies, caused e.g. by attenuation/damping in the inlet/tubing system (Ibrom et al., 2007; Fratini et al., 2012), path averaging, sensor separation, analyser response time, and high- and low-pass filtering; (ii) effects of temperature fluctuations and dilution by water vapour on measured fluctuations in concentrations of CO₂ (Webb–Pearman–Leuning corrections; Webb et al., 1980); and (iii) CO₂ storage below sensor height. Quality assurance and quality control procedures were further developed and agreed upon within CEIP, including statistical tests, non-stationarity, integral turbulence characteristics (Foken et al., 2004) and footprint evaluation (Göckede et al., 2008). Friction velocity (u_∗) threshold filtering was implemented using the moving point test according to Papale et al. (2006) and as described in REddyProc (2019), in order to discard flux data from periods of low turbulence.

Different EC post-processing software was used at the different sites within the project, such that the data were not evaluated in exactly the same way across the CEIP network, but a reasonably good overall agreement was found among the different software, within 5 %–10 % difference for 30 min CO₂ flux values (Mauder et al., 2008; Mammarella et al., 2016). Similarly, for the gap-filling of the 30 min flux time series, during periods of instrument malfunction or unsuitable measurement conditions (low turbulence, insufficient fetch, etc.), and for the partitioning of NEP into GPP and R_eco, a number of alternative algorithms have been developed in the past, based on different sets of principles (Falge et al., 2001; Barr et al., 2004; Reichstein et al., 2005; Lasslop et al., 2010). The gap-filling and partitioning algorithm used by default in this study was the generic online REddyProc (2019) software, implemented also in the European Fluxes Database Cluster. REddyProc was based on (i) Reichstein et al. (2005) for the filling of gaps in the NEE flux data on the basis of information from environmental conditions, (ii) Reichstein et al. (2005) for the night-time-data-based R_eco parameterization (using an Arrhenius-type function of temperature), and (iii) on Lasslop et al. (2010) for the daytime-data-based GPP evaluation (using a rectangular hyperbolic light response curve for NEE and including a temperature sensitivity of respiration and limitation of GPP by vapour pressure deficit).

In this study, for all CEIP flux sites, we have retrieved the fully analysed and validated half-hourly (level-3) and daily to annual (level-4) CO₂ flux (NEP, GPP, R_eco) data as available, initially from the CEIP database and later from the European Fluxes Database Cluster (2012) or from the GHG-Europe portal (GHG-Europe, 2012). For these data, although the evaluation methods were not necessarily harmonized between sites, we hold that the data available in the database were obtained using the best possible, state-of-the-art evaluation methods at the time of retrieval. For the four non-CEIP flux sites, flux evaluation closely followed CEIP protocols; in the case of DB4 the EddyPro (v6.2) software was used, which was based on a synthesis of calculation and correction methods from CEIP and other FLUXNET flux networks around the globe.

The EC-CO₂ flux measurements used in this study mostly spanned the 5-year period of CEIP (2004–2008), except for a dozen sites where measurements continued until 2010, i.e. the end of NEU and of atmospheric N_r sampling. Older EC data (since the mid-late 1990s) were also available at DB5, EN6 and EN13. Data collection started and ended later at DB4, at which both EC-CO₂ flux and DELTA-N_r measurements spanned the 7-year period 2009–2015. Data analyses presented in the paper, based on inter-annual mean CO₂ budgets and mean N_r deposition, assume that five or more years of monitoring yield reasonably robust estimates of long-term fluxes for the different sites and that the small time shift between the CEIP and NEU project periods (2–3 year overlap) does not affect the results significantly. At some sites, such as DB2, long-term NEE measurements showed multi-decadal variations (Pilegaard et al., 2011; Wu et al., 2013); thus it was essential to use the years overlapping with NEU.

2.3.2 Soil CO₂ and CH₄ fluxes

In situ soil CO₂ efflux (SCE) measurements by opaque (static or dynamic) manual chambers were carried out at 25 of the forest sites, with typically weekly to monthly sampling frequency, with fluxes being measured continuously (hourly) by automated chambers at a few sites (e.g. EN2). The SCE is usually considered a proxy for CO₂ production by soil respiration (R_soil), though the two may not be equal as part of the CO₂ production is dissolved into pore water and may reach the atmosphere only later, either on-site or even off-site if dissolved CO₂ (DIC) leaches to groundwater. Annual R_soil data, scaled up from SCE measurements, are available for 19 forest sites and were collected from the CEIP or GHG-Europe databases and/or from various peer-reviewed publications for the different sites (see Table S7). The ratio of heterotrophic respiration (R_het) to R_soil was determined on an annual scale at 16 sites by different techniques (root-exclusion meshes, trenching experiments, radiocarbon or stable isotope tracing, tree girdling; e.g. Subke et al., 2006) (Table S7).

Methane fluxes were measured by chamber methods or eddy covariance at six forest sites and five semi-natural (peatland, wetland) sites (Hendriks et al., 2007; Skiba et al., 2009; Drewer et al., 2010; Shvaleva et al., 2011; Luo et al., 2012; Kowalska et al., 2013; Juszczak and Augustin, 2013) (Table S6). These data were complemented by bioassay measurements of CH₄ emission or uptake (net oxidation) by the laboratory soil cores, as described previously for NO and N₂O estimates (Schauffler et al., 2010).

2.3.3 Dissolved carbon losses

Dissolved inorganic (excluding CO₂ from weathering of carbonate rocks) and organic carbon (DIC/DOC) fluxes were measured at six forest sites (DB1, DB2, EN4, EN8, EN10, EN15), using suction cups for sampling soil water and combined with soil drainage data, or by monitoring water runoff through weirs, as part of CEIP, NEU and other projects (Ilvesniemi et al., 2009; Kindler et al., 2011; Gielen et al., 2011; Verstraeten et al., 2014). Data were also available for peatland at SN7, with DIC, DOC and also dissolved CH₄ concentrations in pore water of the clayey peat, in groundwater from the sand aquifer and in ditch water, as described in Hendriks et al. (2007). For the peatland within SN9, Dinsmore et al. (2010) measured stream concentrations and export of DIC, DOC and particulate organic carbon (POC), and they also estimated stream evasion of CO₂, CH₄ and N₂O in addition to the land-based flux (EC, chamber) measurements in the tower footprint.

2.4 Ecosystem greenhouse gas balance

Net GHG budgets were constructed from inter-annual mean EC-based NEP combined with measured and scaled-up N₂O and CH₄ fluxes wherever available (nine and six sites, respectively), or with bioassay-derived fluxes (most sites) or modelled data (BASFOR, forests/N₂O only), using 100-year global warming potentials (GWPs) of 265 and 28 for N₂O and CH₄, respectively (Fifth Assessment Report, IPCC, 2013). The sign convention for non-CO₂ GHG fluxes and for the net ecosystem GHG balance in this paper adopts an atmospheric warming perspective, i.e. positive fluxes for emissions toward the atmosphere (warming) and negative for uptake by the surface (cooling).

2.5 Ancillary soil, plant and ecosystem measurements

Ancillary data were collected mainly for the purpose of assembling input parameters and calibration datasets for forest ecosystem (BASFOR) modelling (see below). Texture (% clay, % sand, % silt), pH, soil organic carbon (SOC) concentration and C∕N ratios were measured in soils of 35 sites as part of the bioassay experiments described previously but were otherwise also documented in the CEIP database and in papers previously published for the majority of sites (Table S1). For the forest sites, ecosystem data for soil water content (SWC), porosity, saturation water content (Φ_SAT), field capacity (Φ_FC) and wilting point (Φ_WP), as well as for canopy height (H), leaf area index (LAI), diameter at breast height (DBH), basal area (BA), number of trees per unit area or stand density (SD) and thinning events, were obtained from CEIP and other project (e.g. FLUXNET) databases and complemented by various publications (Tables S2–S5). Such was also the case for ecosystem carbon stocks in soil organic matter (CSOM) and in roots (CR), stems (CS), branches (CB), leaves (CL) and litter layers (CLITT), for which the global database assembled by Luyssaert et al. (2007) provided additional data. At sites for which published values of Φ_FC and Φ_WP were not available, default estimates were inferred from soil texture by means of van Genuchten (1980) pedo-transfer functions, using tabulated values from the German soil description handbook (Eckelmann et al., 2005)

Foliar C and N contents (LeafC, LeafN) were measured as part of NEU for EN1, EN2, EN5, EN8, EN10, EN15, EN16, DB2 (Wang et al., 2013), DB4, SN3, SN4, SN8, and SN9 or were otherwise taken from CEIP, GHG-Europe and FLUXNET databases as well as various publications; in total, leaf C∕N measurements were available for 31 sites. By contrast, data were much rarer for C∕N ratios for other compartments of the forest ecosystem, with data available at only 15 sites for litter and at only five sites for roots, stems and branches.

2.6 Modelling of C and N fluxes and pools by the BASFOR ecosystem model

The BASic FORest model (BASFOR) is a deterministic forest ecosystem model that simulates the growth (from planting or natural regeneration) and the biogeochemistry of temperate deciduous and coniferous even-aged stands (van Oijen et al., 2005; Cameron et al., 2013). A description of the model and the fortran code are available in BASFOR (2016). The model was calibrated through a multiple site Bayesian calibration (BC) procedure, applied to three groups of plant functional types (DBF, EN-spruce, EN-pine), based on $\mathrm{C}/\mathrm{N}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ flux and pool data from the CEIP/NEU databases (Cameron et al., 2018). Details on model implementation as part of this study are provided in Flechard et al. (2020).

Briefly, the C, N and water cycles are simulated at a daily time step in interaction with the soil and climate environments and constrained by management (pruning and thinning). Carbon and nitrogen pools are simulated in the different ecosystem compartments (tree stems, branches, leaves and roots, litter layers, and SOM with fast and slow turnover), which are interconnected by internal flows and transformations (e.g. SOM mineralization, nitrogen retranslocation). Carbon, nitrogen and water enter the ecosystem from the atmosphere (photosynthesis, N_r deposition, rainfall). Inorganic nitrogen is taken up from the soil by tree roots; C and N return to the litter and soil pools upon senescence of leaves, branches and roots, and also when trees are pruned or thinned. Losses of C occur through autotrophic (root and shoot) respiration and microbial decomposition into CO₂ of litter and SOM (heterotrophic respiration); losses of N occur through nitrate leaching below the root zone and soil emissions to the atmosphere of NO and N₂O. The water balance is constrained by incoming rainfall, soil water holding capacity and evapotranspiration (ET) simulated by the Penman equation.

3 Results

3.1 Nitrogen inputs and outputs

3.1.1 Nitrogen deposition

Total inorganic N_r deposition ranged from 0.1 to 4.3 g N m⁻² yr⁻¹ across the CEIP/NEU networks (Table 1), with the largest values observed in the Netherlands, northern Belgium, and southern Germany and the lowest levels observed at latitudes > 60^∘ N (Fennoscandia). Nitrogen deposition was dominated by the dry fraction in forests (Fig. 2), with an average contribution to total deposition of 63 % versus 39 % for short semi-natural vegetation. This contribution was even larger (> 2∕3) for high-deposition sites (N_dep > 2 g N m⁻² yr⁻¹). Total N_dep was more strongly correlated to dry deposition across all sites (R²=0.94) than to wet deposition (R²=0.56). Important differences in the ratio of dry to wet deposition are evident across climatic regions, with the share of dry deposition being especially large at Mediterranean sites (e.g. Sanz et al., 2002), where annual rainfall is smaller. However, the share of dry deposition was also large for sites that are located near (large) anthropogenic (industrial, vehicular, agricultural) N_r emission sources. Total N_r deposition was around 25 % smaller on average at short semi-natural vegetation sites compared with forests (Fig. S2), even though the mean total atmospheric N_r concentrations (reduced and oxidized, N-containing gas and aerosol compounds) were quite similar between the two datasets (Flechard et al., 2011). The difference was driven by higher dry deposition rates over forests due to higher aerodynamic roughness and deposition velocities (Fig. S3; see also Schwede et al., 2018). Reduced N_r (NH₃ gas and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in aerosol and rain, collectively NH_x) contributed on average 56 % of total deposition; oxidized N_r (HNO₃+NO₂ gas and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in aerosol and rain, collectively NO_y) was dominant at only six forest sites of the network (EN7, EN10, EN18, EB2, SN3, SN5; Fig. 2).

Figure 2 Total reactive nitrogen deposition (N_dep) and breakdown into inorganic wet and dry, oxidized (NO_y) and reduced (NH_x) deposition estimates at the 31 forest sites (evergreen needleleaf EN1-7 (spruce), EN8-18 (pine), mixed – MF, deciduous broadleaf – DB, evergreen broadleaf – EB) and at nine short semi-natural (SN) vegetation sites of the NitroEurope monitoring network. Data are arithmetic means over the years 2007–2010 of (i) inferential dry deposition estimates by four different models based on in situ atmospheric N_r measurements and of (ii) different wet deposition estimates from precipitation monitoring datasets and from European-scale atmospheric chemistry and transport modelling (EMEP). Error bars indicate standard deviations of the four dry deposition models (red bars) and standard deviations of the different data sources for inorganic N_r wet deposition (blue bars). Wet deposition of water-soluble organic nitrogen (WSON) was measured at a few selected sites and is shown here for comparison with total inorganic N_r deposition.

For comparison, dry deposition, calculated here as the ensemble average of four inferential model estimates based on in situ N_r concentration measurements, was on average more than a factor of 2 larger than the ca. 50 km × 50 km grid square-averaged EMEP model estimate (taken from EMEP, 2013) (see Fig. S2). However, since each EMEP grid square contains variable proportions of different land uses with different deposition velocities, it is more meaningful to compare DELTA-based inferential estimates for each study site with ecosystem-specific EMEP dry deposition rates in the relevant grid squares. In this case, the EMEP dry deposition rates are on average 32 % smaller than the inferential estimates.

By contrast, wet deposition was generally reasonably consistent between the different data sources for inorganic N_r (in situ bulk or wet-only measurement, kriging of monitoring network data, EMEP model output). For the 18 sites where all three sources of data were available, the mean CV of the three estimates was 21 % (range 2 %–56 %, with 15 CV values out of 18 below 30 %), and the mean (±95 % confidence interval) wet deposition estimates across the 18 sites were 0.63±0.14, 0.64±0.15 and 0.68±0.16 g N m⁻² yr⁻¹ for the three methods, respectively (Fig. S2), showing no systematic bias between methods. Wet deposition of organic nitrogen, measured at 16 sites, represented on average 11 % (range 2 %–36 %) of total inorganic + organic wet deposition (Fig. S2), but only 4 % (range 1 %–30 %) of total dry + wet N_r deposition, since total N_dep was dominated by dry deposition at most forest sites.

Figure 3 Comparison of measured and estimated ecosystem inorganic N_r losses and their relationships to total atmospheric N_r deposition (x axis) at the forest sites. NO fluxes (a) and N₂O fluxes (b) were either (i) measured in situ using static or dynamic flux chambers, (ii) scaled up from laboratory-bioassay-derived T–WFPS relationships, or (iii) simulated using the BASFOR ecosystem model (see text for details). DIN leaching (c) was either measured (lysimeter or suction cups) or predicted from the Dise et al. (2009) empirical algorithm. The sum of inorganic N_r losses (d) (DINleach + NO+ N₂O) was computed as the mean of measured values and modelled estimates. In panels (a)–(c), site names are indicated for sites where in situ measurements were available.

3.1.2 Nitrogen losses

Total ecosystem losses of inorganic N_r were computed for the forest sites as the sum of DIN leaching and NO and N₂O emissions (Fig. 3a–d). We assumed that NH₃ emissions by soil and vegetation were negligible due to generally acidic forest soils, as well as low values of the stomatal compensation point (the leaf NH₃ emission potential), respectively (Flechard et al., 2013). Inorganic N_r losses (Fig. 3d) increased sharply with N_r deposition and were largely dominated by DIN leaching at N_dep levels above 2 g N m⁻² yr⁻¹ (Fig. 3c). For these large N_dep levels, the fraction of deposited N_r lost as DIN, NO or N₂O was generally larger than 50 % (Fig. 3f). The inorganic N_r balance (N_r deposition minus NO, N₂O and DIN losses) was probably still positive for most sites (Fig. 3e), although the confidence intervals of the budget term (accounting for uncertainties in all terms including deposition) were very large for the elevated N_r deposition sites. Note that the DIN leaching estimate by BASFOR, shown for comparison on Fig. 3c, was not used in the calculation of total inorganic N losses in Fig. 3d; this is because BASFOR does not simulate N₂ loss by denitrification, and thus part of the soil N surplus that would in reality denitrify is assumed to drain, resulting in an overestimation of the leaching term, though not necessarily of the total N losses.

Emissions of NO estimated from bioassay measurements (Schaufler et al., 2010) and by BASFOR modelling were generally of the same order in forests (average values across all forest sites of 0.22 and 0.21 g N m⁻² yr⁻¹, respectively), but validation by in situ chamber flux data was difficult owing to the limited number of available measurements (only five forest sites, mean value 0.27 g N m⁻² yr⁻¹). Nonetheless, the largest NO emissions by the three methods were all found at N_dep levels above 2 g N m⁻² yr⁻¹ (Fig. 3a). By contrast, N₂O emissions did not show any marked dependence on N_dep and were on average smaller than NO emissions by a factor of 2 to 5, with mean values across all sites of 0.12, 0.08 and 0.04 g N m⁻² yr⁻¹ for bioassay, BASFOR and chamber fluxes, respectively. The mean N₂O fluxes (averaged over the different methods) were larger than mean NO fluxes at only one-third of the forests sites; by contrast, at SN sites N₂O emissions were larger than NO emissions at all but one location. The dominance of NO over N₂O in forests could in principle reflect the generally well aerated conditions of (especially coniferous) forest litter layers on well-drained topsoils, which are more conducive to NO formation by nitrification than N₂O by denitrification (Davidson et al., 2000; Pilegaard et al., 2006). This would be perhaps especially true for the four highest (> 3 g N m⁻² yr⁻¹) N_r deposition sites (EN2, EN8, EN15 and EN16, all coniferous forests) with the highest NO emissions (Fig. 3), which all had sand-dominated (64 %–96 %) soil textures (Table S4). On the other hand, given the acidity of many forest topsoils (Table S4), nitrification could be inhibited, but chemodenitrification could produce significant amounts of NO (Pilegaard, 2013).

For a complete ecosystem net N budget, additional measurements of dissolved organic nitrogen (DON) leaching, as well as dinitrogen (N₂) fluxes (biological fixation and total denitrification), would be required (Fig. 1), but they were not quantified in most cases. A tentative ballpark estimate of the potential magnitude of denitrification N₂ emissions for the DB2 forest site may be calculated by considering the mean N₂∕N₂O ratio of 74 (±0.85 st. err.), which was measured in He–O₂ mixture soil incubation experiments performed on DB2 soil cores (unpublished data). This mean ratio, multiplied by the mean field-measured N₂O emission flux of 0.074 g N₂O-N m⁻² yr⁻¹ (Pilegaard et al., 2006), yields an estimate of the order of 5.5 g N₂-N m⁻² yr⁻¹. There is considerable uncertainty in this number, since the mean N₂∕N₂O ratio was calculated from short-term investigations in the laboratory, which may or may not be representative of the prevailing soil and weather conditions in the field. This uncertainty is reinforced by the low sensitivity of the N₂ detector, which was a factor of 20–80 lower than that of the N₂O detector used in the experiment (Buchen et al., 2019). Another estimate of forest soil denitrification loss obtained through a soil core incubation method was given by Butterbach-Bahl et al. (2002) for the EN2 spruce site, with an annual N₂ emission flux of 0.72 g N₂-N m⁻² yr⁻¹ and a mean N₂/N₂O ratio of 7. The N₂ emissions thus estimated suggest that total denitrification may be a very significant term in the total N budget of forests, possibly of the same order as atmospheric N_r deposition.

Measurements of DON leaching were available at very few sites but proved to be significant. At the pine forest site of EN8, DON leaching was of the order of 0.3 g N m⁻² yr⁻¹, i.e. a factor of 3 lower than DIN losses (Verstraeten et al., 2014). At the beech forest site of DB2, DIN and DON leaching were of the same order (0.07–0.08 g N m⁻² yr⁻¹), but both were very small in comparison to N_dep (2.15 g N m⁻² yr⁻¹), while at the pine forest site of EN10 the leaching/runoff N_r loss was actually dominated by DON (0.012 g N m⁻² yr⁻¹), which was around an order of magnitude larger than DIN leaching (Korhonen et al., 2013) and a factor of 4 smaller than N_dep.

3.2 Net carbon and greenhouse gas balance

3.2.1 Spatial variability of the carbon sink in relation to climate and nitrogen deposition

The ultimate objective of the project was to quantify the response of C sequestration to atmospheric N_r deposition (addressed in Flechard et al., 2020), but this is not straightforward. We follow first in this paper a descriptive approach, whereby variations of C fluxes and other productivity indicators (e.g. leaf area index and N content) are examined graphically as a function of N_dep (Fig. 4). However, this is done with the strong reservation that a simple empirical relationship does not necessarily prove causality, as other confounding and co-varying factors, e.g. climate, soil and age may exist. Figures 4–5 show for example that the large inter-site differences in MAT and MAP at the European scale also need to be considered, beside the variability in N_dep. Note that in assessing the variability of ecosystem carbon sink strength within the network, we use EC-derived NEP (the long-term NEE sum) as a proxy for the net ecosystem carbon balance, because estimates of DIC/DOC leaching, CH₄ emissions and other C loss processes were not systematically measured at all sites.

Figure 4 Overview of inter-annual mean EC-derived C flux estimates (GPP, R_eco and NEP), ecosystem LAI and leaf N content, in relation to total (dry + wet) atmospheric N_r deposition (a–e), and relationship of R_eco to GPP (f), for forests (filled circles, black labels) and short semi-natural vegetation (filled stars, magenta labels). In all plots, the colour scale indicates mean annual temperature (MAT), while the symbol size is proportional to mean annual precipitation (MAP, scale provided in panel a).

Inter-annual mean NEP ranged from a small net source of −70 g C m⁻² yr⁻¹ (EN6, a waterlogged peat-based spruce stand in the southern Russian taiga) to a large net sink of +826 g C m⁻² yr⁻¹ (EN5, upland spruce forest in northern Italy) (Table 1, Fig. 4c); GPP ranged from 377 g C m⁻² yr⁻¹ (SN3, a boreal peatland site with the lowest MAT = −0.6 ^∘C) to 2256 g C m⁻² yr⁻¹ (EN14, a pine stand in Italy, one of the warmest sites with MAT of 14.9  ^∘C and non-limiting rainfall with MAP = 920 mm) (Fig. 4a). Ecosystem respiration peaked at 1767 g C m⁻² yr⁻¹ at EN4 (upland spruce forest in eastern Germany) and was lowest at 345 g C m⁻² yr⁻¹ at SN3 (boreal peatland), the coldest site (Fig. 4b); R_eco was strongly and positively related to GPP (Fig. 4f) (R²=0.62, slope = 0.64). The resulting carbon sequestration efficiency values based on the ratio of observed NEP ∕ GPP (CSE_obs) varied widely among observation sites, ranging from −9 % to 61 %, with an average of 25 %.

The data show a positive correlation between GPP and N_dep in the range 0–2.5 g N m⁻² yr⁻¹ (R²=0.55, p<0.01). By contrast the five sites with N_dep > 2.5 g N m⁻² yr⁻¹ tend to show visually an inverse relationship (Fig. 4a), despite the fact that they lie in comparatively favourable climates. Similar patterns are observed for R_eco and NEP (Fig. 4b–c), but with much larger scatter and lower R² (0.24, p<0.01, and 0.30, p<0.01, respectively, for the N_dep range 0–2.5 g N m⁻² yr⁻¹), with the same apparent decline for higher-deposition sites. However, a closer inspection of Fig. 4a–c reveals a potential cross-correlation with climate (see also Fig. S4): (i) the lower end of the N_dep range, coinciding with the lowest GPP, R_eco and NEP, also coincides with the lowest MAT and MAP (e.g. Finnish sites); and (ii) the sites in the intermediate N_dep range (1.5–2.5 g N m⁻² yr⁻¹), coinciding mostly with the largest observed GPP values (> 1500 g C m⁻² yr⁻¹), were on average 1.8 ^∘C warmer (10.2 vs. 8.4 ^∘C) and 89 mm yr⁻¹ wetter (887 vs. 798 mm) than the sites in the lower N_dep range (0-1.5 g N m⁻² yr⁻¹).

Other proxies of the ecosystem C and N cycles and productivity, such as the LAI (defined as one-sided for broadleaf, or half of the total for needleleaf; Table 1 and Fig. 4d) and the foliar N content (LeafN, Fig. 4e), also showed positive relationships to N_dep (see below for differences between vegetation types). The inter-annual mean value of the annual maximum leaf area index (LAI_max) increased from around 1 to 7 m² m⁻² for N_dep increasing from 0.1 to 4.5 g N m⁻² yr⁻¹, with the lower half of the LAI_max distribution (< 4.5 m² m⁻²) mostly occurring at boreal, Mediterranean, and upland sites and thus under temperature and/or water limitations.

Clearly, therefore, the continental-scale variability in ecosystem–atmosphere CO₂ fluxes was to a large extent controlled by climate, namely by limitations in temperature and water availability. Gross ecosystem productivity was limited, as expected, by low temperatures at high latitudes (or high elevations) and by low rainfall and/or high evaporative demand at Mediterranean, boreal and continental sites. The distribution of the forest monitoring sites in the European climate space, with MAP and MAT on the x and y axes, respectively (Fig. 5a–b), shows that for sites with MAT > 7 ^∘C there was a broad negative correlation between MAT and MAP ($R=-\mathrm{0.49}$, p=0.01); i.e. the warmest sites in southern Europe tend to be the driest and therefore potentially water-limited. Maximum GPP (and also R_eco, not shown) occurred in the mid-climate range, around 9–15  ^∘C MAT and around 700–1000 mm MAP. Similarly, the larger N_dep values (> 2 g N m⁻² yr⁻¹) occurred almost exclusively at sites with MAT in the narrow range of 6–11 ^∘C, and although these large N_dep values were found in a broad MAP range (550–1200 mm), they peaked sharply around 800–900 mm MAP (Figs. 5a; S4). Modelled N_dep values from the EMEP CTM (Fig. 5c–d) show that this is a generic pattern at the European scale.

Figure 5 Distribution of observation-based nitrogen deposition (N_dep) (a) and gross primary productivity (GPP) (b) for the forest sites of this study, within the European climate space represented by mean annual temperature (MAT) and precipitation (MAP). In plot (a) the symbol colour indicates N_dep, while the symbol size is proportional to GPP; in plot (b) the symbol colour indicates GPP, while the symbol size is proportional to N_dep. Plot (c) shows modelled N_dep from the EMEP model over coniferous forests (year 2010), represented in climate space (one data point for each grid square of the EMEP domain containing coniferous forests), also shown as a map (d). The MAT axis can be seen as a proxy for latitude and/or elevation, while the MAP axis expresses to some extent longitude (distance to the ocean) and/or orographic precipitation enhancement.

Ecosystem DIC + DOC losses estimated by Kindler et al. (2011) for four forest sites of this study (DB1, DB2, EN4, EN15) were on average 13±7 g C m⁻² yr⁻¹ (range 3–35 g C m⁻² yr⁻¹), with contributions by DIC to total (DIC + DOC) losses varying between 18 % and 83 %. By contrast, Gielen et al. (2011) estimated DOC leaching losses of 10±2 g C m⁻² yr⁻¹ for the EN8 pine stand on an acidic sandy soil, in which DIC concentrations in soil water were negligibly small. Ilvesniemi et al. (2009) found DOC losses in runoff at EN10 of 0.8 g C m⁻² yr⁻¹, which was negligible compared with NEP. These leaching or runoff losses of DOC and DIC were on average over all forest sites equivalent to a very small mean fraction of 0.6 % of GPP (range 0.1 %–1.9 %) but a more significant fraction of NEP (mean 6 %, range 0.3 %–13 %). At the SN7 peatland site, fluxes of total dissolved carbon (including CH₄) through seepage, infiltration and drainage were relatively small by comparison to NEP and to other peat bogs (17 g C m⁻² yr⁻¹, only 5 % of NEP) (Hendriks et al., 2007); by contrast, at the SN9 peatland site, net stream C export (including DIC, DOC and POC) was on average 29.1 g C m⁻² yr⁻¹ (81 % of which being DOC), equivalent to a mean leached fraction of 37 % of NEP (Dinsmore et al., 2010).

3.2.2 Differences between plant functional types

Forests (F) and short semi-natural (SN) vegetation showed similar relationships with GPP as a function of N_r deposition, increasing with a broadly similar slope at low N_dep values and then levelling off beyond 2g N m⁻² yr⁻¹, except for the fact that GPP was lower by typically 200–500 g C m⁻² yr⁻¹ in SN compared with F sites, for a given N_dep level (Fig. 4). The behaviour was different for NEP, where the slope against N_dep in the range 0–2 g N m⁻² yr⁻¹ was much steeper for F than for SN, which occurred because R_eco values are of the same order for F and SN at a given N_dep level. No systematic difference was observed between the forest PFT, based on the available data, in the apparent relationships of the C fluxes vs. N_dep. However, this may be a result of the small number – and large diversity – of deciduous broadleaf (DB) and evergreen broadleaf (EB) forest sites in the dataset, compared with evergreen needleleaf (EN) sites (Table 1).

The relationship of LeafN to N_dep (Fig. 4e) showed three distinct groups, with the smallest values (0.8 %–1.8 % N in dry weight, DW) for evergreen needleleaf and broadleaf (EN, EB) forests being positively correlated to N_dep in the range 0.5–4.3 g N m⁻² yr⁻¹ (R²=0.71, p< 0.01). Values for short semi-natural (SN) vegetation were found in an intermediate range (1 %–2.7 % N DW), with a steep and significant relationship to N_dep (R²=0.51, p< 0.05). The largest values occurred for deciduous broadleaf (DB) forests (mostly > 2 % N DW), but with little relationship to N_dep (R²=0.18, not significant). Seasonal variations in forest LeafN could reach a factor of 2, as did differences between tree species within the same forest, which may account for some of the scatter observed in Fig. 4e.

Figure 6 BASFOR baseline simulations for all forest sites; model outputs and observation-based values were averaged over the years between the first and last available observations. Note that model simulations include MF and EBF sites, for which the model was not calibrated in Cameron et al. (2018); the two MF runs were made using the parameter table for DBF, while the five EBF runs were made using the parameter table for ENF to allow continued growth throughout the year. H: mean tree height; DBH: mean diameter at breast height; CLBS, CSOM, CR, CLITT: carbon stocks in leaves, branches and stems; in soil organic matter; in roots; and in litter layers, respectively; R²: coefficient of determination; MAE: mean absolute error; NRMSE: root mean square error normalized to the mean.

3.2.3 Carbon fluxes and pools derived from forest ecosystem modelling

In the BASFOR base run (Fig. 6), reasonable overall model performance was achieved for GPP, ecosystem C pools, H, DBH, LAI and LeafN, while more scatter was present for R_eco, NEP and ET. In particular, in apparent contrast to GPP, R_eco stands out as a more challenging variable to model. Predictably, because BASFOR was calibrated using a subset of 22 sites from this dataset (Cameron et al., 2018), the range and mean values of modelled R_eco were close to mean observations by EC across the study sites, but differences between sites were poorly reproduced with much scatter around the 1∕1 line and a low R². The modelled carbon sequestration efficiency (CSE_mod), simulated over the same time period as the flux measurements, was much less variable (range 17 %–31 %, mean 22 %) than observation-based values (CSE_obs) (comparison made for the 22 sites used in model calibration). One possible reason was that BASFOR assumed that autotrophic respiration (R_aut) is a constant fraction of GPP, which may be an oversimplification (Collalti and Prentice, 2019). Also, heterotrophic respiration (R_het) appeared to be a much more variable fraction of R_eco in reality (Table S7) than was predicted by the model, leading to sizable divergence in the overall modelled R_eco. As the direct measurement, NEP was the least uncertain term in EC-derived data, compared with GPP and daytime R_eco, which were inferred from measured (half-hourly) EC NEE by empirical partitioning models. By contrast, in BASFOR, NEP was calculated as the residual between two large numbers (GPP and R_eco) and thus compounds the uncertainties of both component terms. The modelled result for NEP appeared to be an overestimation of net C uptake at low-productivity sites and an underestimation at high-productivity ones (slope < 1). A broadly similar pattern emerged for ET.

Figure 7 Net greenhouse gas (GHG) budgets calculated from a combination of inter-annual mean (around 2005–2010) net ecosystem productivity (NEP) from eddy covariance, and N₂O and CH₄ flux data measured in situ or estimated by extrapolated bioassay data and forest ecosystem BASFOR modelling. Global warming potential values (100-year time horizon) of 265 and 28 were used for N₂O and CH₄, respectively; the sign convention is with respect to the atmosphere: negative for a sink and positive for a source. The data were grouped by ecosystem type (evergreen needleleaf: EN-spruce and EN-pine; MF: mixed forests; DB: deciduous broadleaf; EB: evergreen broadleaf; SN: short semi-natural vegetation); within each group the data were sorted by increasing N_r deposition.

3.2.4 Net ecosystem greenhouse gas budgets

Carbon dioxide largely dominated the net GHG budget at all forest sites, with only three sites where either N₂O or CH₄ GWP-equivalent fluxes were larger than 10 % of NEP in absolute terms (Fig. 7). Most of the forest soils (22 out of 27 sites) investigated in the bioassay experiment behaved as small net sinks for CH₄, with a mean (±SE) net oxidation flux of $-\mathrm{0.14}\pm \mathrm{0.03}$ g C m⁻² yr⁻¹ (range −0.61 to +0.16 g C m⁻² yr⁻¹). The mean CH₄ flux measured by soil chambers at the six forest sites where such measurements were available (EN2, EN6, EN10, EN16, DB2, EB5) was also a net oxidation flux of $-\mathrm{0.32}\pm \mathrm{0.15}$ g C m⁻² yr⁻¹ (range −1.0 to −0.0 g C m⁻² yr⁻¹). For these six sites, there was a significant correlation (R²=0.74, p<0.05) between annual soil CH₄ flux estimates derived from the bioassay experiment and from in situ flux measurements (Fig. S5 in Supplement), with the largest net annual soil CH₄ uptake flux being observed by both methods at the EN10 pine forest site (Skiba et al., 2009). By contrast, at the elevated N_dep sites EN2 and EN16, the net soil CH₄ flux was close to zero, consistent with previous research (e.g. Steudler et al., 1989; Smith et al., 2000) showing that the CH₄ oxidation capacity of forest soils in negatively affected by N_r addition or deposition. In terms of C uptake, soil CH₄ oxidation was negligible compared to CO₂ fluxes, representing on average only 0.1 % of NEP (range 0.0 %–0.4 %). In terms of GWP the CH₄ flux was larger, being equivalent to 0.8 % of NEP (range 0 %–4.5 %), but on average still a factor of 3 smaller than the warming by N₂O emissions equal to 3.9 % of NEP (range 0 %–18.5 %).

In contrast to forests, at semi-natural, short vegetation sites N₂O or CH₄ emissions had a larger impact on the net GHG balance, where most (seven out of nine) sites showed non-CO₂ GHG contributions larger than 10 % of NEP. Three of these seven sites were unfertilized, extensively grazed upland (SN2, SN5, SN6) grasslands (small N₂O sources), while three sites (SN3, SN7, SN8) were CH₄-emitting peatlands or wetlands (EC-CH₄ and chamber flux data from Drewer et al., 2010; Hendriks et al., 2007; Juszczak and Augustin, 2013; Kowalska et al., 2013). At SN3 and SN8, the small to moderate NEP sinks were turned by large CH₄ emissions into net GHG sources (net warming budgets of +127 and +242 g CO₂-C Eq m⁻² yr⁻¹, respectively), though not into actual net C sources (Fig. 7). At SN8, CH₄ emissions generally ranged from 25 to 45 g CH₄-C m⁻² yr⁻¹ but reached 86 g CH₄-C m⁻² yr⁻¹ during a particularly wet year, when the whole area was flooded. At the SN9 peatland site, Dinsmore et al. (2010) calculated that stream GHG evasion – at the scale of the 335 ha peat bog encompassing the flux tower footprint – together with downstream export represented 50–60 g CO₂-Eq m⁻² yr⁻¹ (13–16 g CO₂-C Eq m⁻² yr⁻¹), 96 % of which being degassed CO₂, i.e. in the range 11 %–23 % of the GHG budget from the tower footprint.

4 Discussion

Previous observations of simple empirical relationships found between N deposition and forest productivity have been criticized for, amongst other things, their low number of replications, unreasonably high sensitivities of productivity to N additions, and limitations of the data and simplistic univariate statistical approaches used (Magnani et al., 2007; Högberg, 2007; de Vries et al., 2008; Sutton et al., 2008). Other attempts have subsequently been made to assess impacts of N deposition on forest growth and carbon sequestration, while accounting for other drivers, at more than 350 long-term monitoring plots in Europe (Solberg et al., 2009; Laubhann et al., 2009; De Vries et al., 2008). A special feature of the present study is that it aims to assemble N deposition rates and budgets, together with variables of the carbon cycle, for a wide range of sites across the European continent in more depth and completeness than hitherto attempted, in order to seek more robust empirical evidence for the response of the terrestrial carbon cycle to different regimes of atmospheric N inputs. The quality of the individual datasets is, however, not uniformly high. Some of the data were measured in situ with known uncertainty, while others were simulated, derived from laboratory experiments, and adapted to the field situation using measured time series of soil T and soil moisture, or taken from existing databases and literature. Also, data may not be fully comparable between sites (different methods used) or even fully representative of each site (spatial heterogeneity). In the following sections, we discuss limitations of the measured, empirical and simulated data, in terms of the component C and N fluxes, their budgets and interactions, and the challenges faced when attempting to establish empirical/statistical evidence for possible N effects on carbon sequestration in natural and semi-natural terrestrial ecosystems in Europe.

4.1 Constraining the ecosystem nitrogen balance through combined measurements and modelling

The compilation of N_r flux data (Fig. 3), based on several independent sources for each component term, provides a realistic picture of inorganic N_r inputs and losses; their balance suggests that for forests subjected to large deposition loads (> 2 g N m⁻² yr⁻¹), typically more than half of the incoming N_r is lost to neighbouring environmental compartments such as groundwater and the atmosphere and is thus not available to promote C storage in the forest ecosystem. Since N losses increase – and N retention decreases – exponentially when N_dep exceeds a critical load of approximately 2–2.5 g N m⁻² yr⁻¹ (Fig. 3), it seems unlikely that the C sink strength of semi-natural ecosystems, including forests, increases linearly with N_r deposition, especially not with wet N deposition only. Based on a review of experimental N addition studies (e.g. Högberg et al., 2006; Pregitzer et al., 2008) and monitoring-based field studies along N deposition gradients (e.g. Solberg et al., 2009; Laubhann et al., 2009; Thomas et al., 2010), De Vries et al. (2014) suggested that the C response reaches a plateau near 1.5–2.0 g N m⁻² yr⁻¹ and then starts to decrease. The linear relationship between C sequestration and wet N_r deposition as proposed for example by Magnani et al. (2007) is also challenged by the large contribution of dry N_r deposition and therefore by the poor correlation between total N_dep and wet deposition. We argue that our multiple-constraint approach for the nitrogen balance (measurement–model combination, model ensemble averaging, alternative data sources) provides overall a more robust basis for studying the impact of N_dep on the C cycle, even though uncertainties in individual terms remain significant.

4.1.1 Reducing uncertainty in nitrogen deposition

The uncertainty in dry deposition based on measured N_r concentrations and inferential modelling is likely not smaller than 30 %, due to limitations in process understanding. The difference between ecosystem-specific EMEP values and the mean inferential estimates (Fig. S2) reflects discrepancies and uncertainties in the four dry deposition schemes used (Flechard et al., 2011); the mean coefficient of variation (CV = σ∕μ) between the four inferential model estimates was 36 %, i.e. larger than the difference between ecosystem-specific EMEP values and the mean inferential estimates. Other sources of discrepancy between the two methods include the use of measured vs. modelled meteorology to drive the deposition models, and site-specific vs. generic values of canopy height and leaf area index, as discussed in Flechard et al. (2011).

The uncertainty in total N_r deposition is probably of the same order since even wet deposition can be deceptively difficult to measure (Dämmgen et al., 2005), and organic N, especially water-soluble organic N, may be significant but challenging to quantify (Cape et al., 2012) and generally ignored in the literature. WSON appears to be a generally small fraction of total (wet + dry) N_dep at most sites except at remote locations in Fennoscandia (EN10, SN3), where WSON deposition could represent up to 20 %–30 % of total N_dep. Also, potential double-counting due to dry deposition to the bulk deposition collectors (e.g. Thimonier et al., 2018) was not considered in this study, although on the basis of the comparison to other data sources (Fig. S2), bulk samplers did not appear to significantly overestimate wet deposition.

Despite these uncertainties, measuring gas-phase and aerosol N_r concentrations locally should provide a better estimate of total ecosystem N_r inputs than the outputs of a large-scale chemical transport model. In addition, the partitioning of wet vs. dry deposition, reduced vs. oxidized N, and canopy absorption vs. soil deposition should also be improved, all of which are useful in interpreting ecosystem N-cycling processes. In particular, for ammonia, with its high spatial variability on a local scale, the inferential modelling approach based on local measurements is likely to provide more realistic deposition estimates than a coarse-resolution chemical transport model (Flechard et al., 2013; Thimonier et al., 2018). In addition to low-cost methods for N_r concentrations, more actual micrometeorological N_r flux measurements are needed to further process understanding and better constrain surface exchange models over many ecosystems (Fowler et al., 2009). For example, ammonia flux measurements at DB2 have revealed unexpected features such as net NH₃ emissions from the forest in summer and autumn, in particular in response to leaf fall (Hansen et al., 2013, 2017). DB2 is likely not a net NH₃ source at the annual scale, but short-term emission pulses, which are not represented in most dry deposition models (Flechard et al., 2011), could significantly offset total N_r deposition.

An improved knowledge of N_r exchange patterns over CO₂ flux monitoring sites, either through inferential modelling or direct flux measurements, is also essential to quantify the fraction of deposited N_r that is absorbed by the canopy, reaching more or less directly the seat of photosynthesis in leaves, thus favouring a higher nitrogen use efficiency (NUE) (Nair et al., 2016; Wortman et al., 2012; Gaige et al., 2007). Canopy nitrogen retention occurs via several processes, including gaseous uptake by stomatal diffusion, a well-documented process (Monteith and Unsworth, 1990), but also through cuticular diffusion and stomatal penetration by aqueous solutions, with surface-deposited and dissolved gases and particles acting as direct leaf nutrients (Burkhardt, 2010; Burkhardt et al., 2012). By contrast, the N_r fraction initially deposited to soil (as simulated by the majority of fertilization tracer experiments, e.g. Nadelhoffer et al., 1999) is subject to various losses via nitrification, denitrification and microbial uptake, before being eventually taken up by roots and moving upwards in xylem flow. The more advanced, emerging multi-layer canopy exchange models for atmospheric pollutants (N_r species but also O₃, SO₂, etc.) can now partition dry deposition into stomatal, non-stomatal and soil pathways with increasing detail (Zhou et al., 2017; Simpson and Tuovinen, 2014; Flechard et al., 2013), thanks to improved understanding and parameterizations of surface and air column interactions and of photosynthesis-driven stomatal conductance (Büker et al., 2007; Grote et al., 2014). However, particular attention must be paid to measurement quality for an improved deposition accuracy, because such models are still very much dependent on local atmospheric concentration data for all main N_r forms (gas and aerosol, reduced and oxidized, mineral and organic).

4.1.2 Uncertainty in ecosystem nitrogen losses and net balance

The comparison of DIN leaching values by different methods shows that the Dise et al. (2009) algorithm performs reasonably well for low to moderate N_r deposition levels but underestimates DIN losses for some of the highest (> 4 g N m⁻² yr⁻¹) deposition sites. This observation was also made by Dise et al. (2009) themselves, who argued that their simple relationships involving external forcings (N_dep) and internal factors (soil N status) are adequate “for early to intermediate stages of nitrogen saturation” but may fail at sites where historical, chronically enhanced N_r deposition has so strongly impacted forest ecosystems that N leaching has become dependent also on stochastic factors such as insect defoliation or a drought period followed by re-wetting of the soil. As was the case for field-measured NO emissions (Fig. 3a), the four highest DIN leaching fluxes (0.9–3.2 g N m⁻² yr⁻¹) occurred in the four highest N_dep forests growing on well-drained acidic sandy soils. In addition, it is noteworthy that the two sites with the largest N_dep and DIN leaching rates (EN15, EN16) were dominated by pine or Douglas fir (Table S2). These species have been shown in a common garden experiment (Legout et al., 2016) to cause larger nitrification, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and acidification rates (as well as larger losses of calcium, magnesium and aluminium), compared with other tree species such as beech or oak. This is consistent with deciduous trees being known to take up and store more nitrogen per unit biomass in stems and branches than coniferous trees (Jacobsen et al., 2003). Typical stem N content values, proposed for N uptake calculations in the Convention on Long-range Transboundary Air Pollution (CLRTAP) manual for critical loads mapping, are 1 and 1.5 g N kg⁻¹ dry matter for conifers and deciduous trees, respectively, for steady-state conditions (CLRTAP, 2017). Tree species traits may therefore, in our study, have exacerbated an existing DIN leaching predisposition resulting from edaphic factors and pollution climate. At the lower end of the N_dep range, the dataset is consistent with previous studies, which have shown that DIN leaching is unlikely to occur in forests where N_dep < 1 g N m⁻² yr⁻¹ (de Vries et al., 2009), although under these conditions there may still be significant N losses as NO and N₂O (Fig. 3).

The best empirical fit for the relationship of the sum DIN + NO + N₂O to N_dep was slightly non-linear (Fig. 3d) and may indicate that at the upper end of the N_dep range, above 4 g N m⁻² yr⁻¹, the sum of inorganic N_r losses might approach or even exceed the estimated atmospheric deposition, which corresponds to one of the several existing definitions of ecosystem N saturation (see below). Whether these ecosystems turn into net N sources depends on the relative magnitudes of the missing terms: N₂ fixation (likely small in temperate compared with tropical forests; Vitousek et al., 2002), N₂ losses from denitrification (possibly the largest of the unknown terms at forest sites that are frequently waterlogged), N₂O losses from the litter layers of the forest floor, DON leaching, and also incoming organic nitrogen in precipitation (WSON) as well as dry deposition of organic N_r species, not quantified here (Fig. 1). The presumably small, and unaccounted for, N inputs via N₂ fixation and organic N_r deposition are at least partly compensated for by denitrification N₂ losses and DON leaching losses. Moreover DON leaching typically responds much less strongly than DIN leaching to N inputs (Siemens and Kaupenjohann, 2002). Under these assumptions, the inorganic N_r budget calculated from Fig. 3 may provide a reasonable proxy for the overall ecosystem N balance. In this case, N outputs by gaseous and dissolved losses represent on average across all forest sites 43 % of N inputs. More important than the average N loss for judging N_r deposition effects on C sequestration is the large range of losses from 6 % to 85 %, with on average a 27 % loss (range 6 %–54 %) for N_dep < 1 g N m⁻² yr⁻¹, 45 % loss (12 %–78 %) for intermediate N_dep levels, and 65 % loss (35 %–85 %) for N_dep > 3 g N m⁻² yr⁻¹. However, if the very few available data or estimates for DON leaching and especially denitrification N₂ fluxes are correct and may be extrapolated to other sites, they may often outweigh the inputs through organic N_r deposition and biological N₂ fixation, and thus the inorganic N_r budget (Fig. 3) may underestimate the overall N losses.

4.2 Drivers and uncertainties of the carbon and GHG balance

4.2.1 Variability of carbon sequestration efficiency

The CSE ratio (= NEP ∕ GPP) calculated over the CEIP/NEU project observation periods provides an indicator of the fraction of accumulated carbon in the ecosystem relative to gross CO₂ uptake by photosynthesis. This is a useful metric to compare carbon cycling in different terrestrial ecosystems, and it is directly related to climate effects and other drivers such as site fertility (Vicca et al., 2012) and management (Campioli et al., 2015). By contrast, quantifying the accumulated carbon in terrestrial ecosystems requires much longer observations (one or several decades) to ensure statistical significance of a small change over a large C stock, particularly when soils are considered. This is often impractical but also of limited use, because N deposition rates are unlikely to be constant over such long periods.

Over the time frame of the CEIP/NEU projects, observation-based CSE_obs values were much more variable than their modelled counterparts. Negative CSE_obs values (EN6, EN11) imply a net carbon source and may be explained by a number of factors, including soil carbon loss, lateral DOC ∕ DIC water flow from adjacent ecosystems, tree mortality, low fertility, poor ecosystem health, a recently planted forest or other disturbances with long-lasting consequences on the C budget. At EN6 the main reasons may be a large SOC concentration, leading to large R_eco values, and a relatively old age of the forest, responsible for a small GPP.

However, the large discrepancy between observation-based and modelled CSE estimates may not be entirely caused by the model's inability to reproduce all fine patterns of GPP and especially R_eco across all ecosystems (Fig. 6). Some of the largest CSE_obs values may be less ecologically plausible and might result from methodological biases and/or incorrect interpretation of the EC measurements, in terms of their representativeness for the ecosystem considered. Multi-annual values of GPP and R_eco derived from EC flux data are not measurements sensu stricto; they compound problems in EC measurements, post-processing of high-frequency data, gap-filling and partitioning. Some partitioning algorithms (Barr et al., 2004; Reichstein et al., 2005) evaluate GPP as the difference between measured daytime NEE and an estimate of daytime R_eco that is based on an empirical model of night-time R_eco measurements. In this case, any problem with night-time and thus with estimated daytime R_eco would directly impact GPP in the same way (Vickers et al., 2009): GPP and R_eco would both be underestimated, or both overestimated, in absolute terms and by the same absolute magnitude, thereby impacting the annual or long-term NEP ∕ GPP (CSE_obs) ratio.

In this study, however, the use of the daytime-data-based partitioning method by Lasslop et al. (2010), within the REddyProc algorithm embedded in the European Fluxes Database Cluster, was intended to ensure the independence of GPP and R_eco estimates, since R_eco was estimated from the intercept of the Michaelis–Menten light response curve fitted to daytime-measured NEE. This partitioning procedure should avoid the propagation into the GPP estimate of potential errors in night-time R_eco data, although it still assumes similar dependencies of day- and night-time respiration on environmental factors, which is debatable from a biological standpoint (e.g. Kok, 1949; Wehr et al., 2016; Wohlfahrt and Galvagno, 2017). From a micrometeorological perspective, the night-time flux can be underestimated due to low-turbulence conditions and the transport of CO₂ by horizontal and/or vertical advection, as well as the decoupling of soil-level and understorey fluxes from the turbulent fluxes measured above the canopy (Feigenwinter et al., 2008; Etzold et al., 2010; Montagnani et al., 2010; Paul-Limoges et al., 2017). Further, in principle, the u_∗ threshold filtering (Gu et al., 2005; Papale et al., 2006), carried out to discard low-turbulence flux data at the start of the gap-filling and partitioning algorithm (REddyProc, 2019) should alleviate the issue of night-time R_eco underestimation, which affects annual R_eco and CSE_obs even if the error does not propagate into GPP in the Lasslop et al. (2010) method. However, the choice of the value for the u_∗ threshold can be critical if advection-affected flux values are to be discarded, especially for sites and datasets where the independence of the gap-filled annual NEP value from the u_∗ threshold value cannot be demonstrated. Advective flux contributions remain a largely unresolved issue, as Aubinet et al. (2010) conclude that “direct advection measurements do not help to solve the night-time CO₂ closure problem”. Others (e.g. Kutsch and Kolari, 2015) have commented on the need to assign appropriate uncertainties when dealing with CSE and C balances derived from EC flux towers, which only measure turbulent fluxes and CO₂ storage change in the air column underneath the sensor but not the other terms of the conservation equation of a scalar in the atmospheric boundary layer (see Eq. 1 in Aubinet et al., 2000).

Despite all these precautions, at sloping or complex terrain sites where advection can be important, it cannot be excluded that the Lasslop et al. (2010) daytime-data-based approach may still underestimate R_eco (and overestimate CSE_obs) if advection is not accounted for explicitly. This is because the R_eco estimate based on the intercept of the light response curve for the measured NEE (at PAR = 0) is strongly influenced by measurements made around sunrise and sunset, when a clear impact of advection on the light response curve ordinate has been observed, as shown at the EN5 subalpine site by Montagnani et al. (2009) (see their Fig. 13).

It is important to note that advection may also be a problem at flat lowland sites if there is strong spatial land surface heterogeneity, e.g. differences in albedo or in Bowen ratio, a gradient in tree species, a nearby lake, and a gradient in water availability. Conversely, there may also be sites where EC underestimates CSE_obs for similar reasons, albeit in the opposite direction, for example additional CO₂ being advected into the ecosystem and then released by turbulent diffusion to the atmosphere within the tower footprint. Another possibility is that basal R_eco, measured at dawn or dusk over a different (larger) footprint, is lower than during the day. Flux partitioning may again in this case underestimate R_eco during the warmer daytime hours, and therefore also underestimate GPP, resulting in overestimated NEP ∕ GPP (CSE_obs) ratios.

Figure 8 Variability of observation-based and modelled carbon sequestration efficiency (CSE, defined as the NEP ∕ GPP ratio), as a function of (a) site elevation above mean sea level (m), and (b) MAP: mean annual precipitation (mm). Site labels are provided for observations only.

Given this uncertainty, the fact that most of the forest stands with CSE_obs values larger than 40 % (EN1, EN5, DB6, MF2) were located at elevations above 700 m a.s.l. (Table 1 and Fig. 8a), i.e. in hilly or mountainous areas with topographically more complex terrain than typically encountered at lowland sites, may be coincidental, or partly a consequence of advection or decoupling issues (Paul-Limoges et al., 2017). In such conditions, consistency cross-checks involving additional flux, advection, soil and biometric measurements, even ecosystem modelling, provide useful reference points to assess the plausibility of EC-derived C budgets and to better constrain the problem. At the EN5 site, the annual total tree biomass C increment based on biometric measurements was on average 218 g C m⁻² yr⁻¹ over the period 2010–2017 (Leonardo Montagnani, unpublished data), i.e. 26 % of the reported mean EC-derived NEP value of 826 g C m⁻² yr⁻¹ for the CEIP–NEU period, and it seems unlikely that the increase in soil carbon and fine root stocks could account for the large difference. By contrast, the DB6 site was a fertile and managed beech forest, with a significantly higher efficiency conversion of photosynthates into biomass compared to less fertile and unmanaged sites (Vicca et al., 2012; Campioli et al., 2015). The long-term annual total NPP at the site was 780 g C m⁻² yr⁻¹ over the period 1992–2007, with a significant part allocated below ground (Alberti et al., 2015), while heterotrophic respiration estimated at the site using either bomb carbon (Harrison et al., 2000) or mineralization rates (Persson et al., 2000) was around 200 g C m⁻² yr⁻¹, resulting in similar NEP estimates by EC flux measurements versus biometric data combined with process studies.

At the MF2 site, Etzold et al. (2011) calculated inter-annual mean EC-derived NEP, GPP and R_eco values (for the same 2005–2009 period used in this study) of 415, 1830 and 1383 g C m⁻² yr⁻¹, respectively, using customized gap-filling and partitioning algorithms and thus providing alternative estimates to those from the REddyProc algorithm within the European Fluxes Database Cluster (Table 1). Values of R_eco and NEP were 82 % larger and 40 % lower, respectively, in Etzold et al. (2011) compared with the default database values that do not explicitly correct for advection. However, the Etzold et al. (2011) mean EC-derived NEP was much closer to NEP values calculated from the net annual increment in the woody and non-woody biomass and soil C storage using four different biometric and modelling methods (range 307–514 g C m⁻² yr⁻¹, mean 421 g C m⁻² yr⁻¹). The CSE_obs value derived from Etzold et al. (2011) was 23 % and comparable to the value of 25 % that can be calculated from the decoupling-corrected EC budget computed by Paul-Limoges et al. (2017) for the same site for the years 2014–2015, in which the decoupling correction to account for undetected below-canopy fluxes doubled R_eco and reduced NEP from 758 to 327 g C m⁻² yr⁻¹. These alternative CSE_obs estimates were thus much lower than the default CSE_obs value of 48 % but fully consistent with model predictions (Fig. 8).

The four upland sites, EN1, EN5, DB6 and MF2, were also among the wettest, with MAP > 1000 mm (Fig. 8b) in principle promoting larger leaching and runoff. The overall distribution of CSE_obs as a function of MAP (Fig. 8b) shows an apparent increase in CSE_obs with precipitation, though with large scatter, which would be consistent with a reduction in EC-tower-based R_eco through an increase in the dissolved leached fraction. At sites where significant leaching occurs, R_eco determined from the atmospheric flux is no longer a reliable indicator of total C losses by respiration since the dissolved and then leached fraction of R_soil is not captured by the flux tower (Gielen et al., 2011), which implies that CSE_obs is overestimated. As observed in the case of GPP, such apparent correlations of CSE_obs to single factors like elevation or MAP may not be (entirely) causal, potentially concealing underlying cross-correlations (such as large but unmeasured advection components occurring at the same sites where MAP is largest). The data by Kindler et al. (2011) and Gielen et al. (2011) do suggest that the overestimation of C sequestration (as estimated by EC-derived NEP), caused by not accounting for dissolved C leaching, was likely smaller than 10 % for forests (7 % of NEP on average), but all five sites they investigated had MAP < 1000 mm and only one (EN4) was an upland site (785 m).

To summarize a set of unresolved issues, the largest CSE_obs values (> 45 %) are likely to result from a combination of ecological factors and methodological biases, but they occurred at sites in the mid-range for N_dep (1.2–2.2 g N m⁻² yr⁻¹) and thus did not introduce confounding trends in the overall C–N relationships we seek to establish across the whole N_dep spectrum in this study.

4.2.2 Forest net greenhouse gas balance dominated by carbon

Based on the available data, the net GHG balance of the 31 forests investigated was generally not significantly affected by N₂O or CH₄ (Fig. 7), with the caveat that these fluxes were not actually measured in situ everywhere or with the same intensity and duration as CO₂. Thus, the uncertainty in non-CO₂ GHG fluxes is much larger (possibly > 100 %) than for multi-annual EC-based CO₂ datasets, where a typical uncertainty is of the order of 10 %–30 % (Loescher et al., 2006). Nonetheless, the N₂O and CH₄ emissions observed by different methods in forest soils were typically 2 orders of magnitude smaller than the CO₂ sink (in GWP equivalents), which means that the quality of CO₂ estimates dominates the overall uncertainty in our forest GHG budgets. Note that such results cannot be extended to waterlogged, organic soils of temperate and boreal zones, where CH₄ emissions can be large (Morison et al., 2012), nor can they be extended to the tropics, especially in degraded forests (Pearson et al., 2017). Also, N₂O fluxes can be highly episodic, with emission events linked to for example freeze–thaw cycles (Risk et al., 2013; Medinets et al., 2017), and such episodes would have been missed by the bioassay approach and in some cases by discontinuous (manual) chamber measurements.

By contrast, for the short semi-natural vegetation sites of our study, NEP was on average a factor of 2.7 smaller than in forests but only a factor of 1.5 smaller for GPP, which implies that total C losses were much larger in proportion to gross assimilation, especially non-respiratory, non-CO₂ losses (i.e. a much lower CSE). Large wetland CH₄ emissions and dissolved DIC/DOC fluxes were much more likely to offset or even determine the C and GHG balance (Fig. 7; Kindler et al., 2011). In these systems, studying the impact of N_r deposition on C sequestration requires much more robust estimates of the gaseous and dissolved budgets for all components and over the long term, since the estimation of NECB requires in addition to EC CO₂ the knowledge of non-atmospheric, non-CO₂ fluxes (Fig. 1). Technological developments in the field of (routine) EC measurements for N₂O and CH₄ (e.g. Nemitz et al., 2018) are likely to reduce uncertainties in net GHG budgets in the foreseeable future, but DIC/DOC losses in wetlands probably represent a bigger challenge.

It should however be remembered that such short-term GHG budgets, based on a few years of flux data and GWP multipliers for a 100-year time horizon, do not actually reflect the long-term climate impact of northern mires, which may be thousands of years old and, despite their CH₄ emissions, typically have an overall climate-cooling effect. As shown by Frolking et al. (2006), pristine mires typically start cooling the climate some hundreds of years after their formation, with the exact timing of course depending on the magnitude of the CH₄ and CO₂ fluxes; thus the history of the site should be accounted for when dealing with ecosystem radiative forcing assessments. For the SN3 site, Drewer et al. (2010) actually used a 500-year time horizon GWP (instead of the usual 100-year) for CH₄, reducing the GHG source strength of the site by a factor of 4 to 10, depending on the year considered.

Figure 9 Relationships of leaf (a) and topsoil (b) C∕N ratios with atmospheric nitrogen deposition (N_dep), and to each other (c), in different ecosystem types (DBF: deciduous broadleaf forests, MF: mixed forests, ENF: evergreen needleleaf forests, EBF: evergreen broadleaf forests, SN: short semi-natural vegetation).

4.3 Challenges in understanding the coupling of carbon and nitrogen budgets

4.3.1 Tangled effects of nitrogen deposition and climate on ecosystem productivity

The analysis of N_dep variability and spatial patterns at the scale of the monitoring network, as well as the European scale (Fig. 5), showed that the impact of N_r deposition on ecosystem C sequestration cannot be considered independently of climate in the regional context of this study. Nitrogen deposition patterns at the European scale result from the continent-wide geographical distribution of population, human, industrial and agricultural activities, and of precursor emissions, combined with mesoscale patterns of meteorology-driven atmospheric circulation and chemistry. Through the interplay of these factors, the elevated N_dep levels in this study happened to co-occur geographically with temperate climatic zones of central and western Europe (Fig. 5c–d) that are the most conducive to vegetation growth at the continental scale. This means adequate water supply as precipitation, reasonably low summertime evaporative demand, mild winters and temperate summers, and long growing seasons. In other words, there are many gaps in the multidimensional variable space, which is incompletely explored by the available dataset. Thus, any regression analysis that would correlate NEP and other C fluxes with N_dep, without simultaneously accounting for climate, would be flawed, as Sutton et al. (2008) concluded from their reanalysis of the data used by Magnani et al. (2007). A dC∕dN slope calculated directly from a (linear or non-linear) mono-factorial regression analysis of GPP or NEP vs. N_dep would misleadingly attribute the whole C flux variability to N_dep while ignoring climate effects (Fleischer et al., 2013). In addition, a range of other potential explanatory variables such as soil type, especially the water holding capacity (Φ_FC−Φ_WP), soil fertility (Vicca et al., 2012; Legout et al., 2014), tree species and stand age (Besnard et al., 2018), are potentially needed to explain the observed variability. In order to account for, and untangle, the multiple inter-relationships, we chose a mechanistic model (BASFOR) based approach, described in Flechard et al. (2020), whereby most of the known interactions of plant, soil, climate, age and species are encoded and parameterized to the best of our current knowledge. Given the limited size and very large diversity of the dataset, such an approach appears to be preferable to regression-based statistical analyses, since a simple pattern to explain the coupling of carbon and nitrogen budgets with the available data and knowledge is unlikely.

4.3.2 Evidence of nitrogen saturation from various indicators

Various definitions of nitrogen saturation have been proposed (Aber, 1992; De Schrijver et al., 2008; Binkley and Högberg, 2016), including (i) the absence of a growth response in the case of further N addition (dC∕dN = 0), (ii) the onset of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and/or gaseous emissions, and (iii) the equivalence of N inputs and N losses. The underlying concept of a dC∕dN response is that the C and N cycles are closely coupled through stoichiometric ratios in the different parts of the ecosystem, with very different C∕N ratios in soil organic matter, roots, leaves, tree branches and stems (de Vries et al., 2009; Zechmeister-Boltenstern et al., 2015). A difference in dC∕dN response could, for example, be expected between forests, where carbon is stored in both woody and root biomass (C∕N ratio 300–500) and below ground in SOM (C∕N ratio 30–40), versus short semi-natural vegetation, where most of the stock is in SOM, and thus with a much lower overall ecosystem C∕N ratio. This would be consistent with the observations in Fig. 4, where the apparent increase in NEP with increasing N_dep is smaller in short semi-natural vegetation than in forests. But the theoretical stoichiometric approach becomes more uncertain in the event of N saturation, as the C and N cycles have become much less tightly coupled than in pristine, N-limited environments, and thus defining a dose-response relationship requires a precise quantification of all C and N inputs and losses, not just productivity and N_r deposition.

Another possible indicator of N saturation in the present dataset may be provided by the comparison of the relationships of C∕N ratios of foliage and topsoil (5 cm) to atmospheric N_r deposition (Fig. 9a–b). Since leaf N content was not only dependent on N_r deposition but also on the ecosystem type (Fig. 4e), C∕N ratios are shown separately for the different vegetation classes in Fig. 9. There was a clear negative correlation of leaf C∕N ratio to N_dep for coniferous forests (ENF, spruce and pine pooled: exponential fit R²=0.86, p< 0.01) and a similar but not significant trend for SN (linear R²=0.29) (Fig. 9a); for the other ecosystems (DBF, MF, EBF) there were not enough data to derive trends. In topsoils (Fig. 9b), there was also a broad downward trend of C∕N ratios with increasing N_dep within the ENF and SN classes but only for N_dep up to 2.5 g N m⁻² yr⁻¹. Again, as for GPP and NEP, the relationship is highly non-linear as the four ENF sites above this N_dep threshold break the trend observed in the lower N_dep sites, and the overall best fit is quadratic (R²=0.49, p< 0.01) with an inflexion point around this threshold. While the relationship of foliar C∕N ratio to N_dep was almost linear for ENF (a consequence of the linear trend in ENF leaf N content, Fig. 4e), the non-linear behaviour of the topsoil C∕N ratio and its stabilization or increase for N_dep > 2.5 g N m⁻² yr⁻¹ indicate a possible threshold for saturation. Atmospheric nitrogen was therefore apparently efficiently taken up by vegetation when reaching the leaves; but after leaf fall, and following litter decomposition and incorporation into the topsoil, there appeared to be a limit to the amount of nitrogen that can be stabilized into soil organic matter of the ENF sites. However, forest soil organic N stocks are very large (in the range 200–700 g N m⁻² at the sites we investigated), and therefore changes in C∕N ratios in response to atmospheric N_r deposition must be very slow. The soil C∕N ratio at a given time reflects centuries of land use as well as a more recent history of multi-decadal changes in N_r deposition (Flechard et al., 2020). This complicates the interpretation of the downward trends observed from instantaneous snapshots of soil and foliar C∕N ratios versus N_dep since the ecosystems cannot be considered to be in steady state, neither for N_dep nor for growth or productivity. There was a positive correlation across all vegetation types between topsoil and foliar C∕N ratios (Fig. 9c; R²=0.19, p< 0.05), but this was mostly driven by differences between plant functional types (no significant correlation within each PFT).

Following definition (ii) of N saturation given above, the sum of inorganic N_r losses, heavily dominated by DIN leaching at the upper end of the N_dep range in our datasets (Fig. 3), may indicate various stages of N saturation in all forests with N_dep > 1–1.5 g N m⁻² yr⁻¹. A threshold for a more advanced saturation stage could be placed at 2–2.5 g N m⁻² yr⁻¹, where inorganic N_r losses are consistently larger than 50 % of N_dep. Such numbers are entirely consistent with the leaching risk classification of European forests proposed by Dise and Wright (1995), with low leaching risk at N_dep < 1 g N m⁻² yr⁻¹, intermediate risk at N_dep in the range 1–2.5 g N m⁻² yr⁻¹ and high risk at N_dep > 2.5 g N m⁻² yr⁻¹. The results are also in line with the review by De Vries et al. (2014); based on literature results of dC∕dN responses derived from stoichiometric scaling, meta-analysis of N addition experiments, and field observations of both growth changes and N_r deposition, accounting for other drivers, the data showed beneficial N_r deposition effects up to 2–3 g N m⁻² yr⁻¹ and adverse effects at higher levels. A lower N_dep threshold of 1 g N m⁻² yr⁻¹ had also been suggested by de Vries et al. (2007), but this was using throughfall deposition, which generally underestimates total deposition through canopy retention processes (Thimonier et al., 2018). It must be stressed, however, that the definition of an all-purpose, generic N_dep threshold for N saturation may be misleading, or at least qualified with an uncertainty, since some tree species (Douglas fir, pine, spruce), grown on the same soil and under the same climate and N_dep regime, may result in significantly higher ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching rates than others (Legout et al., 2016). This also means that the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching flux is not necessarily a good proxy of the severity of N saturation, though this depends on which of the several definitions of N saturation is considered.

The upper threshold of 2–2.5 g N m⁻² yr⁻¹ happens to coincide with the levelling off of GPP, R_eco and NEP, as well as the further reduction in C fluxes at higher N_dep levels (Fig. 4a-c). Whether this should be interpreted as a negative impact of advanced N saturation on soil processes and plant functioning and, hence, C sequestration potential is not straightforward (Binkley and Högberg, 2016). If the parallel effects of climate, soil fertility, other nutrient limitations, tree species traits, age and planting density are overlooked in a simplistic, first-order interpretation, the dataset hints at an optimum N_dep level around 2 g N m⁻² yr⁻¹, beyond which no further benefits (in carbon terms) could be gained from further atmospheric N_r additions, which would be consistent with the 2–2.5 g N m⁻² yr⁻¹ N_dep threshold derived by Etzold et al. (2014) for Swiss forests. The high soil N_r losses observed in these ecosystems growing under relatively favourable climates would then suggest that whatever fertilization effect N_r deposition may have at low to moderate deposition rates (< 2 g N m⁻² yr⁻¹) is unlikely to be sustained at high-deposition levels, especially on acidic sandy soils. However, the very limited number of affected sites with N_dep > 3 g N m⁻² yr⁻¹ leaves too few degrees of freedom to make the argument statistically compelling. More importantly, a knowledge of all other limitations to growth (climate, soil, fertility, nutrients, age structure) would be required to confirm the hypothesis. Additional measurement- and model-based investigations to untangle the N_dep effect on C sequestration (the dC∕dN term) are presented in Flechard et al. (2020), drawing from the results, fluxes and budgets presented here.

5 Conclusions

We provided estimates of carbon, nitrogen and greenhouse gas budgets for 40 flux tower sites over European forests and semi-natural vegetation, compiled from a large variability of state-of-the-art methods that can be applied in such a network approach. The CO₂ budgets from well-established EC methods were the least uncertain, followed by GHG budgets of forests and then the CH₄ and DIC/DOC fluxes of wetlands; uncertainty levels were likely highest in the net N budgets, especially at the elevated N_r deposition sites where ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching was almost of the same order as N_dep. The uncertainty was still compounded by the lack of some data on biological N₂ fixation; N₂ loss by denitrification; and organic N_r in rainwater, in dry deposition and in soil leaching, but some of these unknown terms would compensate mutually to some extent. Nevertheless, the low-cost network to monitor atmospheric gas-phase and aerosol N_r contributed to substantially reducing the large uncertainty in total N_dep rates at individual sites (compared with gridded outputs of a regional chemical transport model). This was because dry deposition almost systematically heavily dominates over wet deposition in forests, except at very remote sites (away from sources of atmospheric pollution), and directly measured N_r concentrations reduced the uncertainty in dry deposition fluxes.

The greenhouse gas balances of the 31 forest sites included in this study were almost entirely determined by the CO₂ budgets, with small to negligible contributions by N₂O and CH₄. The GHG balance of nine extensively managed and upland grasslands, moorlands and wetlands was much more dependent on CH₄ and N₂O fluxes. Ecosystem productivity (GPP, NEP) data across Europe showed an apparent increase with atmospheric N_dep, though only up to 2.5 g N m⁻² yr⁻¹, while the larger N_dep rates also happen to coincide geographically with regions of Europe where climate is optimal for tree growth (neither too cold nor too dry). The data thus underpinned a strong covariation of N_r deposition with variables like elevation and climate, and they indicated that the ecosystem response of carbon sequestration to nitrogen deposition cannot be calculated simply and directly from the observed apparent dNEP ∕ dN_dep using bivariate statistics. Other co-varying influences such as climate, soil, fertility, nutrient availability, forest age and ecophysiological processes should be analysed alongside so the nitrogen deposition effect can be isolated.

The site-specific analysis of C and N fluxes and budgets across a large geographical and climatic gradient support the concept of a non-linear response of C sequestration to N deposition. Large nitrogen losses (especially nitrate) from forests suggest that up to one-third of the sites investigated can be classified as in the early to advanced stages of N saturation. At the sites with the largest N_r deposition rates (> 2.5 g N m⁻² yr⁻¹), a stagnation or reduction in forest productivity, compared to mid-range deposition sites, was observed. Beyond the conclusion that the apparent C response to increased N_r deposition was non-linear, we do not have enough data to test the hypothesis that the reduction in productivity and C sequestration is linked to N-saturation-induced ecological impacts on soil and ecosystem functioning, rather than just the confounding effects of variability in meteorological and other drivers. Further efforts are required to disentangle N_dep effects and climatic as well as pedological effects on C sequestration at the continental scale.