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 pe-riodical 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 ﬂux towers across Europe to local measurement-based estimates of dry and wet N r deposition from a dedicated col-located 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 ﬂux measurements of NO, N 2 O and CH 4 ﬂuxes; soil NO − 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 (BAS-FOR) modelling.

Abstract. 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 2 O and CH 4 fluxes; soil NO − 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 (BAS-FOR) modelling.
Multi-year averages of net ecosystem productivity (NEP) in forests ranged from −70 to 826 g C m −2 yr −1 at total wet + dry inorganic N r deposition rates (N dep ) of 0.3 to 4.3 g N m −2 yr −1 and from −4 to 361 g C m −2 yr −1 at N dep rates of 0.1 to 3.1 g N m −2 yr −1 in short semi-natural vegetation (moorlands, wetlands and unfertilized extensively managed grasslands). The GHG budgets of the forests were strongly dominated by CO 2 exchange, while CH 4 and N 2 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 2 losses by denitrification. Nitrogen losses in the form of NO, N 2 O and especially NO − 3 were on average 27 % (range 6 %-54 %) of N dep at sites with N dep < 1 g N m −2 yr −1 versus 65 % (range 35 %-85 %) for N dep > 3 g N m −2 yr −1 . 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 −2 yr −1 , 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 −2 yr −1 ), where in-organic 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 .

Introduction
The global terrestrial net sink for atmospheric carbon dioxide (CO 2 ) is approximately 1.7 Pg C yr −1 , i.e. roughly onefifth of global CO 2 -C emissions by fossil fuel combustion and industry (9.4 ± 0.5 Pg C yr −1 ). This corresponds to the land-based carbon (C) uptake of 3.2 ± 0.8 Pg C yr −1 minus emissions from deforestation and other land-use changes of 1.5 ± 0.7 Pg C yr −1 . The ocean sink is of the same order (2.4 ± 0.5 Pg C yr −1 ), while twice as much CO 2 -C (4.7 ± 0.02 Pg C yr −1 ) is added yearly to the atmosphere (Le Quéré et al., 2018). Data from atmospheric CO 2 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 2 sink is currently occurring in forests, and especially in boreal and temperate forests of the Northern Hemisphere 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 2 sink may be stronger -and the northern hemispheric land CO 2 sink weaker -than was generally believed. At the European scale, Schulze et al. (2010) calculated that the net biome productivity (NBP, the mean longterm 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 2 sinks have been attributed to a combination of factors including afforestation of abandoned land and formerly cut forests, reduced forest harvest, CO 2 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(Magnani et al., , 2008Hö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  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 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) 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 , 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. 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 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 (NO − 3 ) leaching (Dise et al., 2009), nitric oxide (NO) and nitrous oxide (N 2 O) emissions (Pilegaard et al., 2006); ammonia (NH 3 ) bidirectional exchange (Hansen et al., 2013); and emissions of di-nitrogen (N 2 ) from total denitrification (Butterbach-Bahl et al., 2002) (Fig. 1). The implication is that the carbon response to N dep would be nonlinear, 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 (NH + 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).
Carbon losses through dissolved organic carbon (DOC) and biogenic dissolved inorganic carbon (DIC) leaching can also be significant, especially for wetlands  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 4 ) may further represent a significant loss from organic soils (Hendriks et al., 2007), while CH 4 oxidation, which is often observed in well-aerated soils and can be suppressed by N r addition, especially NH + 4 (Steudler et al., 1989), may affect the net greenhouse gas (GHG) budget. Nitrogen-deposition-induced N 2 O emissions from the forest floor (Pilegaard et al., 2006;Liu and Greaver, 2009), or from denitrification triggered by deposited NO − 3 in peatland , 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 2 respiration (Treseder, 2008); a reduction of basal respiration without significant decline in total microbial biomass, following N addition to incubated peat cores ; and added NO − 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 4 emissions from wetlands and peatlands 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 2 flux data. The net ecosystem carbon balance (NECB) includes in addition other C loss fluxes such as DIC/DOC, CH 4 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 2 O. The total N budget includes in addition organic nitrogen deposition (WSON) and leaching (DON), as well as N 2 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). , since in anaerobic conditions and in the presence of NO − 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 4 emissions (Nykänen et al., 2002), as the aerenchyma of tracheophytes provides a direct diffusion path to the atmosphere for soil-produced CH 4 , 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 NO − 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.

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.
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 Mediter-ranean 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 semiarid) 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 2 flux measurement stations, most sites (36) were part of the CEIP CO 2 flux network. A further three CO 2 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.

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.

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 3 . 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.
Atmospheric inorganic N r concentrations, available from the NEU (2013) database, were measured monthly for 2-4 years in the gas phase (NH 3 , HNO 3 , HONO) and in the aerosol phase (NH + 4 , NO − 3 ), using DEnuder for Long-Term Atmospheric sampling (DELTA) systems (Sutton et al., 2001;Tang et al., 2009). Concentrations of nitrogen dioxide (NO 2 ), 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)(2008)(2009)(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 Ncontaining gas and aerosol species, whereby measured am-     Kindler et al. (2011). 10 Kowalska et al. (2013). 11 Legout et al. (2016). 12 Luo et al. (2012). 13 Pilegaard et al. (2006). 14 REddyProc (2019). 15 Schaufler et al. (2010). 16 Simpson et al. (2012). 17 Tang et al. (2009). 18 van Oijen et al. (2005). 19 See Table S7. bient N r concentrations were multiplied by a vegetation-, meteorology-and chemical-species-dependent deposition velocity (V d ) (Flechard et al., , 2013Bertolini et al., 2016;Thimonier et al., 2018). In the case of NH 3 , 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 (NH + 4 , NO − 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 measurementbased 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)(2009)(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 (NH + 4 and NO − 3 ) as well as SO 2− 4 , Cl − , PO 3− 4 , base cations (Mg 2+ , Ca 2+ , 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 2 concentrations, discussed above), simulated by the European-scale EMEP CTM (Simpson et al., 2006a(Simpson et al., , b, 2012 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 2 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.

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 2 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 2 O at seven forest sites and four seminatural sites, as well as at five forest sites for NO (Table S6). Manual static chamber N 2 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 2 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 2 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 2 GHG gas flux measurements were unavailable at many sites. These soil N 2 O and NO (and also CH 4 ) 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 waterfilled 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 2 , N 2 O, NO and CH 4 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 = NH + 4 + NO − 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 −2 yr −1 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.

Ecosystem-atmosphere CO 2 exchange
Half-hourly rates of net ecosystem-atmosphere CO 2 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 4 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 2 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, mediumterm 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 2 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 2 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 2 concentration and the vertical component of wind speed (w), and blockaveraging 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 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 2 (Webb-Pearman-Leuning corrections; Webb et al., 1980); and (iii) CO 2 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 2 flux values 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 REd-dyProc (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-databased 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 2 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 2 flux measurements used in this study mostly spanned the 5-year period of CEIP (2004)(2005)(2006)(2007)(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 2 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 2 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 Wu et al., 2013); thus it was essential to use the years overlapping with NEU.

Soil CO 2 and CH 4 fluxes
In situ soil CO 2 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 2 production by soil respiration (R soil ), though the two may not be equal as part of the CO 2 production is dissolved into pore water and may reach the atmosphere only later, either on-site or even off-site if dissolved CO 2 (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).

Dissolved carbon losses
Dissolved inorganic (excluding CO 2 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 4 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 2 , CH 4 and N 2 O in addition to the land-based flux (EC, chamber) measurements in the tower footprint.

Ecosystem greenhouse gas balance
Net GHG budgets were constructed from inter-annual mean EC-based NEP combined with measured and scaled-up N 2 O and CH 4 fluxes wherever available (nine and six sites, respectively), or with bioassay-derived fluxes (most sites) or modelled data (BASFOR, forests/N 2 O only), using 100year global warming potentials (GWPs) of 265 and 28 for N 2 O and CH 4 , respectively (Fifth Assessment Report, IPCC, 2013). The sign convention for non-CO 2 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).

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.

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 C/N/H 2 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 en-vironments 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 2 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 2 O. The water balance is constrained by incoming rainfall, soil water holding capacity and evapotranspiration (ET) simulated by the Penman equation.

Nitrogen deposition
Total inorganic N r deposition ranged from 0.1 to 4.3 g N m −2 yr −1 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 −2 yr −1 ). Total N dep was more strongly correlated to dry deposition across all sites (R 2 = 0.94) than to wet deposition (R 2 = 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 . 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 3 gas and NH + average 56 % of total deposition; oxidized N r (HNO 3 + NO 2 gas and NO − 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).
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 −2 yr −1 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.

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 2 O emissions ( Fig. 3a-d). We assumed that NH 3 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 3 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 −2 yr −1 (Fig. 3c). For these large N dep levels, the fraction of deposited N r lost as DIN, NO or N 2 O was generally larger than 50 % (Fig. 3f). The inorganic N r balance (N r deposition minus NO, N 2 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 2 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 −2 yr −1 , 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 −2 yr −1 ). Nonetheless, the largest NO emissions by the three methods were all found at N dep levels above 2 g N m −2 yr −1 (Fig. 3a). By contrast, N 2 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 −2 yr −1 for bioassay, BASFOR and chamber fluxes, respectively. The mean N 2 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 2 O emissions were larger than NO emissions at all but one location. The dominance of NO over N 2 O in forests could in principle reflect the generally well aerated conditions of (especially coniferous) forest litter layers on welldrained topsoils, which are more conducive to NO formation by nitrification than N 2 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 −2 yr −1 ) N r deposition sites (EN2, EN8, EN15 and EN16, all coniferous forests) with the highest NO emissions (Fig. 3), which all had sanddominated (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 2 ) 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 2 emissions for the DB2 forest site may be calculated by considering the mean N 2 /N 2 O ratio of 74 (±0.85 st. err.), which was measured in He-O 2 mixture soil incubation experiments performed on DB2 soil cores (unpublished data). This mean ratio, multiplied by the mean field-measured N 2 O emission flux of 0.074 g N 2 O-N m −2 yr −1 (Pilegaard et al., 2006), yields an estimate of the order of 5.5 g N 2 -N m −2 yr −1 . There is considerable uncertainty in this number, since the mean N 2 /N 2 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 2 detector, which was a factor of 20-80 lower than that of the N 2 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 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.
al. (2002) for the EN2 spruce site, with an annual N 2 emission flux of 0.72 g N 2 -N m −2 yr −1 and a mean N 2 /N 2 O ratio of 7. The N 2 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 −2 yr −1 , 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 −2 yr −1 ), but both were very small in comparison to N dep (2.15 g N m −2 yr −1 ), while at the pine forest site of EN10 the leaching/runoff N r loss was actually dominated by DON (0.012 g N m −2 yr −1 ), which was around an order of magnitude larger than DIN leaching ) and a factor of 4 smaller than N dep .

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 4 emissions and other C loss processes were not systematically measured at all sites. Inter-annual mean NEP ranged from a small net source of −70 g C m −2 yr −1 (EN6, a waterlogged peat-based spruce stand in the southern Russian taiga) to a large net sink of +826 g C m −2 yr −1 (EN5, upland spruce forest in northern Italy) (Table 1, Fig. 4c); GPP ranged from 377 g C m −2 yr −1 (SN3, a boreal peatland site with the lowest MAT = −0.6 • C) to 2256 g C m −2 yr −1 (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 −2 yr −1 at EN4 (upland spruce forest in eastern Germany) and was lowest at 345 g C m −2 yr −1 at SN3 (boreal peatland), the coldest site (Fig. 4b); R eco was strongly and positively related to GPP (Fig. 4f)  The data show a positive correlation between GPP and N dep in the range 0-2.5 g N m −2 yr −1 (R 2 = 0.55, p<0.01). By contrast the five sites with N dep > 2.5 g N m −2 yr −1 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 2 (0.24, p<0.01, and 0.30, p<0.01, respectively, for the N dep range 0-2.5 g N m −2 yr −1 ), with the same apparent decline for higherdeposition 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 −2 yr −1 ), coinciding mostly with the largest observed GPP values (> 1500 g C m −2 yr −1 ), were on average 1.8 • C warmer (10.2 vs. 8.4 • C) and 89 mm yr −1 wetter (887 vs. 798 mm) than the sites in the lower N dep range (0-1.5 g N m −2 yr −1 ).
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 2 m −2 for N dep increasing from 0.1 to 4.5 g N m −2 yr −1 , with the lower half of the LAI max distribution (< 4.5 m 2 m −2 ) 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 2 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 = −0.49, p = 0.01); i.e. the warmest sites in south- 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).
ern 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 −2 yr −1 ) 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.
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 −2 yr −1 (range 3-35 g C m −2 yr −1 ), 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 −2 yr −1 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 −2 yr −1 , 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 4 ) through seepage, infiltration and drainage were relatively small by comparison to NEP and to other peat bogs (17 g C m −2 yr −1 , 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 −2 yr −1 (81 % of which being DOC), equivalent to a mean leached fraction of 37 % of NEP .

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 −2 yr −1 , except for the fact that GPP was lower by typically 200-500 g C m −2 yr −1 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 −2 yr −1 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 −2 yr −1 (R 2 = 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 2 = 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 2 = 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.

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 2 . 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.

Net ecosystem greenhouse gas budgets
Carbon dioxide largely dominated the net GHG budget at all forest sites, with only three sites where either N 2 O or CH 4 GWP-equivalent fluxes were larger than 10 % of NEP in absolute terms (Fig. 7). Most of the forest soils (22 out 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.
of 27 sites) investigated in the bioassay experiment behaved as small net sinks for CH 4 , with a mean (±SE) net oxidation flux of −0.14 ± 0.03 g C m −2 yr −1 (range −0.61 to +0.16 g C m −2 yr −1 ). The mean CH 4 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 −0.32±0.15 g C m −2 yr −1 (range −1.0 to −0.0 g C m −2 yr −1 ). For these six sites, there was a significant correlation (R 2 = 0.74, p<0.05) between annual soil CH 4 flux estimates derived from the bioassay experiment and from in situ flux measurements (Fig. S5 in Supplement), with the largest net annual soil CH 4 uptake flux being ob-served by both methods at the EN10 pine forest site . By contrast, at the elevated N dep sites EN2 and EN16, the net soil CH 4 flux was close to zero, consistent with previous research (e.g. Steudler et al., 1989;Smith et al., 2000) showing that the CH 4 oxidation capacity of forest soils in negatively affected by N r addition or deposition. In terms of C uptake, soil CH 4 oxidation was negligible compared to CO 2 fluxes, representing on average only 0.1 % of NEP (range 0.0 %-0.4 %). In terms of GWP the CH 4 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 2 O emissions equal to 3.9 % of NEP (range 0 %-18.5 %). 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 2 : coefficient of determination; MAE: mean absolute error; NRMSE: root mean square error normalized to the mean. In contrast to forests, at semi-natural, short vegetation sites N 2 O or CH 4 emissions had a larger impact on the net GHG balance, where most (seven out of nine) sites showed non-CO 2 GHG contributions larger than 10 % of NEP. Three of these seven sites were unfertilized, extensively grazed upland (SN2, SN5, SN6) grasslands (small N 2 O sources), while three sites (SN3, SN7, SN8) were CH 4 -emitting peatlands or wetlands (EC-CH 4 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 4 emissions into net GHG sources (net warming budgets of +127 and +242 g CO 2 -C Eq m −2 yr −1 , respectively), though not into actual net C sources (Fig. 7). At SN8, CH 4 emissions generally ranged from 25 to 45 g CH 4 -C m −2 yr −1 but reached 86 g CH 4 -C m −2 yr −1 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 2 -Eq m −2 yr −1 (13-16 g CO 2 -C Eq m −2 yr −1 ), 96 % of which being degassed CO 2 , i.e. in the range 11 %-23 % of the GHG budget from the tower footprint.

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 productiv-ity 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 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.

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 −2 yr −1 ), 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 −2 yr −1 (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 −2 yr −1 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.

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 ; the mean coefficient of variation (CV = σ/µ) between the four inferential model estimates was 36 %, i.e. larger than the difference between ecosystemspecific 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 chal-lenging 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 largescale 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 . For example, ammonia flux measurements at DB2 have revealed unexpected features such as net NH 3 emissions from the forest in summer and autumn, in particular in response to leaf fall (Hansen et al., 2013(Hansen et al., , 2017. DB2 is likely not a net NH 3 source at the annual scale, but short-term emission pulses, which are not represented in most dry deposition models , could significantly offset total N r deposition. An improved knowledge of N r exchange patterns over CO 2 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 3 , SO 2 , 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).

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 −2 yr −1 ) 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 −2 yr −1 ) 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, NO − 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 −1 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 −2 yr −1 , although under these conditions there may still be significant N losses as NO and N 2 O (Fig. 3).
The best empirical fit for the relationship of the sum DIN + NO + N 2 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 −2 yr −1 , 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 2 fixation (likely small in temperate compared with tropical forests; Vitousek et al., 2002), N 2 losses from denitrification (possibly the largest of the unknown terms at forest sites that are frequently waterlogged), N 2 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 2 fixation and organic N r deposition are at least partly compensated for by denitrification N 2 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 −2 yr −1 , 45 % loss (12 %-78 %) for intermediate N dep levels, and 65 % loss (35 %-85 %) for N dep > 3 g N m −2 yr −1 . However, if the very few available data or estimates for DON leaching and especially denitrification N 2 fluxes are correct and may be extrapolated to other sites, they may often outweigh the inputs through organic N r deposition and biological N 2 fixation, and thus the inorganic N r budget (Fig. 3) may underestimate the overall N losses.

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 2 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 observationbased 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. Multiannual 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 lowturbulence conditions and the transport of CO 2 by horizontal and/or vertical advection, as well as the decoupling of soillevel 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 lowturbulence 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 advectionaffected 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 2 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 2 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 2 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.
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 in- volving 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 −2 yr −1 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 −2 yr −1 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 −2 yr −1 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 −2 yr −1 , 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 interannual 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 −2 yr −1 , respectively, using customized gapfilling 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 −2 yr −1 , mean 421 g C m −2 yr −1 ). 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 −2 yr −1 . 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 −2 yr −1 ) 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.

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 2 O or CH 4 (Fig. 7), with the caveat that these fluxes were not actually measured in situ everywhere or with the same intensity and duration as CO 2 . Thus, the uncertainty in non-CO 2 GHG fluxes is much larger (possibly > 100 %) than for multi-annual EC-based CO 2 datasets, where a typical uncertainty is of the order of 10 %-30 % (Loescher et al., 2006). Nonetheless, the N 2 O and CH 4 emissions observed by different methods in forest soils were typically 2 orders of magnitude smaller than the CO 2 sink (in GWP equivalents), which means that the quality of CO 2 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 4 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 2 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 2 losses (i.e. a much lower CSE). Large wetland CH 4 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 2 the knowledge of non-atmospheric, non-CO 2 fluxes (Fig. 1). Technological developments in the field of (routine) EC measurements for N 2 O and CH 4 (e.g.  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 4 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 4 and CO 2 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 4 , reducing the GHG source strength of the site by a factor of 4 to 10, depending on the year considered. 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.

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 NO − 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 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 conif-1610 C. R. Flechard et al.: Carbon-nitrogen interactions in European ecosystems -Part 1 erous forests (ENF, spruce and pine pooled: exponential fit R 2 = 0.86, p< 0.01) and a similar but not significant trend for SN (linear R 2 = 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 −2 yr −1 . 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 2 = 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 nonlinear behaviour of the topsoil C/N ratio and its stabilization or increase for N dep > 2.5 g N m −2 yr −1 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 −2 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 2 = 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 −2 yr −1 . A threshold for a more advanced saturation stage could be placed at 2-2.5 g N m −2 yr −1 , 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 −2 yr −1 , intermediate risk at N dep in the range 1-2.5 g N m −2 yr −1 and high risk at N dep > 2.5 g N m −2 yr −1 . 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, metaanalysis 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 −2 yr −1 and adverse effects at higher levels. A lower N dep threshold of 1 g N m −2 yr −1 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 NO − 3 leaching rates than others (Legout et al., 2016). This also means that the NO − 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 −2 yr −1 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 −2 yr −1 , 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 −2 yr −1 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 −2 yr −1 ) 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 −2 yr −1 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 measurementand 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.

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 2 budgets from well-established EC methods were the least uncertain, followed by GHG budgets of forests and then the CH 4 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 NO − 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 2 fixation; N 2 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 2 budgets, with small to negligible contributions by N 2 O and CH 4 . The GHG balance of nine extensively managed and upland grasslands, moorlands and wetlands was much more dependent on CH 4 and N 2 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 −2 yr −1 , 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 −2 yr −1 ), 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-saturationinduced 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.
Code and data availability. The data used in this study are publicly available from online databases and from the literature as described in the "Materials and methods" section.
The codes of models and other software used in this study are publicly available online as described in the "Materials and methods" section.
Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/bg-17-1583-2020-supplement.  Review statement. This paper was edited by Sönke Zaehle and reviewed by A. Àvila and one anonymous referee.