2.2.1 General description of the BASFOR ecosystem model

The BASic FORest (BASFOR) model is a process-based, deterministic forest ecosystem model, which simulates the growth and biogeochemistry (C, N and water cycles) of temperate deciduous and coniferous stands at a daily time step (van Oijen et al., 2005; Cameron et al., 2013, 2018). Model code and documentation are available on GitHub (BASFOR, 2016). Interactions with the atmospheric and soil environments are simulated in some detail, including the role of management (thinning or pruning). BASFOR is a one-dimensional model, i.e. no horizontal heterogeneity of the forest is captured, and BASFOR does not simulate some variables which are important in forest production, such as wood quality or pests and diseases.

Nine state variables for the trees describe (i) C pools – leaves, branches and stems; roots; and reserves (CL, CB and CS or collectively CLBS; CR; CRES; kg C m⁻²); (ii) N pool in leaves (NL; kg N m⁻²); and (iii) stand density (SD, trees m⁻²), tree phenology (only for deciduous trees), accumulated chill days (d) and accumulated thermal time (Tsum; ^∘C d). Seven state variables for the soil can be divided into three categories, according to the three biogeochemical cycles being simulated: (i) C pools in litter layers of the forest floor (CLITT), soil organic matter (SOM) with fast turnover (CSOMF) and SOM with slow turnover (CSOMS) (kg C m⁻²); (ii) N pools as for C but also including mineral N (NLITT, NSOMF, NSOMS, NMIN; kg N m⁻²); and (iii) the water pool: amount of water to the depth of soil explored by the roots (WA; kg H₂O m⁻² = mm) (see Table 1).

Table 1BASFOR model state variables, inputs and outputs, and other acronyms used in the study.

Download XLSX

Carbon enters the system via photosynthesis, calculated as the product of photosynthetically active radiation (PAR) absorption by the plant canopy and light use efficiency (LUE). The leaf and branch pools are subject to senescence, causing carbon flows to litter. Roots are also subject to senescence, causing a flow to fast-decomposing soil organic matter. Litter carbon decomposes to fast-decomposing soil organic matter plus respiration. Fast-decomposing soil organic matter decomposes to slow-decomposing soil organic matter plus respiration. Finally, the slow organic carbon pool decomposes very slowly to CO₂. Nitrogen enters the system in mineral form through atmospheric deposition. Nitrogen leaves the system through leaching and through emission of N₂O and NO from the soil to the atmosphere. N₂ losses from denitrification and biological N₂ fixation are not simulated. Dissolved inorganic nitrogen (DIN) is taken up by the trees from the soil, and nitrogen returns to the soil with senescence of leaves, branches and roots, and also when trees are pruned or thinned. Part of the N from senescing leaves is reused for growth. The availability of mineral nitrogen is a Michaelis–Menten function of the mineral nitrogen pool and is proportional to root biomass. The model does not include a dissolved organic nitrogen (DON) pool and therefore does not account for the possible uptake of bioavailable DON forms (e.g. amino acids, peptides) by trees. Transformations between the four soil nitrogen pools are similar to those of the carbon pools, with mineral nitrogen as the loss term. Water is added to the soil by precipitation and lost through transpiration, evaporation and drainage. Evaporation and transpiration are calculated using the Penman equation, as functions of the radiation intercepted by the soil and vegetation layer, and atmospheric temperature, humidity and wind speed. Drainage of ground water results from water infiltration exceeding field capacity of the soil.

In BASFOR, the C and N cycles are coupled in both trees and soil. The model assumes that new growth of any organ proceeds with a prescribed N∕C ratio, which is species-specific but generally higher for leaves and roots than for stems and branches. If the nitrogen demand for growth cannot be met by supply from the soil, some of the foliar nitrogen is recycled until leaves approach a minimum N∕C ratio when leaf senescence will be accelerated. The calculation of foliar senescence accounts for a vertical profile of nitrogen content, such that the lowest leaves have the lowest N∕C ratio and senesce first. Nitrogen deficiency, as measured by foliar nitrogen content, not only increases leaf senescence, but also decreases GPP and shifts allocation from leaves to roots. Given that foliar N content is variable in BASFOR, the litter that is produced from leaf fall also has a variable N∕C ratio. When the litter decomposes and is transformed, the N∕C ratio of the new soil organic matter will therefore vary too in response to the ratio in the litter. Except for woody plant parts, the C and N couplings in BASFOR vegetation and soil are based on the same generic ecophysiological assumptions as those explained in detail for the grassland model BASGRA (Höglind et al., 2020).

The major inputs to the model are daily time series of weather variables (global radiation, air temperature, precipitation, wind speed and relative humidity). The last two of these are used in the calculation of potential rates of evaporation and transpiration. Soil properties, such as parameters of water retention (field capacity, wilting point, soil depth), are provided as constants. Further, the model requires time series indicating at which days the stand was thinned or pruned. The model outputs include, amongst others, the state variable for trees and soil as well as evapotranspiration (ET), groundwater recharge, canopy height (H), leaf area index (LAI), diameter at breast height (DBH), GPP, R_eco and R_soil, NEP, N mineralization, N leaching, and NO and N₂O emissions (Table 1).

2.2.2 Model implementation and calibration

BASFOR simulations of forest growth and C, N and H₂O fluxes were made for all CEIP–NEU forest sites from planting (spanning the interval 1860–2002) until the end of the NEU project (2011). At a few sites, natural regeneration occurred, but for modelling purposes a planting date was assigned based on the age of the trees. Meteorological data measured at each site over several years since the establishment of the flux towers (typically 5–10 years) were replicated backwards in time in order to generate a time series of model inputs for the whole period since planting. Assumptions were made that inter-annual meteorological variability was sufficiently covered in the span of available measurements and that the impact of climate change since planting was small and could be neglected.

The atmospheric CO₂ mixing ratio was provided as an exponential function of calendar year, fitted to Mauna Loa data since the beginning of records in 1958 (NOAA, 2014) and extrapolated backwards to around 1860–1900 for the oldest forests included in this study. The global CO₂ mixing ratio driving the model thus increased from around 290 ppm in 1900 to 315 ppm in 1958 to 390 ppm in 2010 (Fig. 1). Similarly, atmospheric N_r deposition was a key input to the model and was forced to vary over the lifetimes of the planted forests; N_dep was assumed to rise from pan-European levels well below 0.5 g N m⁻² yr⁻¹ at the turn of the 20th century, sharply increasing after World War II to reach an all-time peak around 1980, and decreasing subsequently from peak values by about one-third until 2005–2010, at which point the NEU N_dep estimates were obtained. We assumed that all sites of the European network followed the same relative time course of N_dep over the course of the 20th century, taken from van Oijen et al. (2008) but scaled for each site using the NEU N_dep estimates (Fig. S1 in the Supplement).

Forest management was included as an input to the model in the form of a prescribed time course of stand density and thinning from planting to the present date. Tree density was known at all sites around the time of the CEIP–NEU projects (Table S2 in Flechard et al., 2020), but information on thinning history since planting (dates and fractions removed) was much sparser. A record of the last thinning event was available at only one-third of all sites, and knowledge of the initial (planting) density and a reasonably complete record of all thinning events were available at only a few sites. For the purposes of BASFOR modelling, we attempted to recreate a plausible density and thinning history over the lifetime of the stands. The guiding principle was that after the age of 20 years one could expect a decadal thinning of the order of 20 %, following Cameron et al. (2013), while the initial reduction was 40 % during the first 20 years. In the absence of an actual record of planting density (observed range: 1400–15 000 trees ha⁻¹), a default initial value of 4500 trees ha⁻¹ was assumed (for around two-thirds of the sites). The general principles of this default scheme were then applied to fit the available density and thinning data for each site, preserving all actual data in the time series while filling in the gaps by plausible interpolation. The density time courses thus obtained, underlying all subsequent model runs, are shown in Fig. S2.

The model was calibrated through a multiple site Bayesian calibration (BC) procedure, applied to three groups of plant functional types (PFTs), based on $\mathrm{C}/\mathrm{N}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ flux and pool data from the CEIP–NEU databases (see Cameron et al., 2018). A total of 22 sites were calibrated, including deciduous broadleaf forests (DB1–6) and evergreen needleleaf forest ENF-spruce (EN1–7) and ENF-pine (EN8–18). The model parameters were calibrated generically within each PFT group; i.e. they were not optimized or adjusted individually for each observation site. In the companion paper (Flechard et al., 2020), baseline BASFOR runs were produced for all 31 forest sites of the network, including also those stands for which the model was not calibrated, such as Mediterranean evergreen broadleaf (EB1 through EB5) and mixed deciduous–coniferous (MF1, MF2), to test the predictive capacity of the model beyond its calibration range (see Fig. 6 in Flechard et al., 2020). However, for the analyses and scenarios presented hereafter, these seven uncalibrated sites were removed from the dataset, as were two additional sites: EN9 and EN12 – EN9 because this agrosilvopastoral ecosystem called “dehesa” has a very low tree density (70 trees ha⁻¹; Tables S1–S2 in the Supplement to Flechard et al., 2020) and is otherwise essentially dry grassland for much of the surface area, which BASFOR cannot simulate; EN12 because this was a very young plantation at the time of the measurements, also with a very large fraction of measured NEP from non-woody biomass. All the conclusions from BASFOR meta-modelling are drawn from the remaining 22 deciduous, pine and spruce stands (sites highlighted in Table S1).

2.2.3 Modelling time frames

In the companion paper (Flechard et al., 2020), C and N budgets were estimated primarily on the basis of ecosystem measurements and for the time horizon of the CEIP and NEU projects (2004–2010). In this paper, BASFOR simulations of the C and N budgets for the 22 forest sites were considered both (i) over the most recent 5-year period (around the time of CEIP–NEU), which did not include any thinning event and started at least 3 years after the last thinning event (referred to hereafter as “5-year”); and (ii) over the whole time span since forest establishment, referred to here as “lifetime”, which ranged from 30 to 190 years across the network and reflected the age of the stand at the time of the CEIP–NEU projects. Note that the term “lifetime” in this context was not used to represent the expected age of senescence or harvest.

On the one hand, the short-term (5-year) simulations were made to evaluate cases where no disturbance by management impacted fluxes and pools over a recent period, regardless the age of the stands at the time of the C and N flux measurements (ca. 2000–2010). On the other hand, the lifetime simulations represent the time-integrated flux and pool history since planting, which reflects the long-term C sequestration (NECB) potential, controlled by the cumulative impact of management (thinning), increasing atmospheric CO₂ mixing ratio, and changing N_r deposition over the last few decades. Thinning modifies the canopy structure and therefore light, water and nutrient availability for the trees, and reduces the LAI momentarily, and in theory the leftover additional organic residues (branches and leaves) could increase heterotrophic respiration and affect the NEP. However, the impact of the disturbance on NEP and R_eco is expected to be small and short-lived (Granier et al., 2008), and a 3-year wait after the last thinning event appears to be reasonable for the modelling. The 5-year data should in theory reflect the C∕N flux observations, although there were a few recorded thinning events during the CEIP–NEU measurement period, and the thinning sequences used as inputs to the model were reconstructed and thus not necessarily accurate (Fig. S2).

2.2.4 Modelled carbon sequestration efficiency (CSE) and nitrogen uptake efficiency (NUPE)

For both C and N, we define modelled indicators of ecosystem retention efficiency relative to a potential input level, which corresponds to the total C or N supply, calculated over both 5-year (no thinning) and lifetime horizons to contrast short-term and long-term patterns. For C sequestration, the relevant terms are the temporal changes in carbon stocks in leaves, branches and stems (CLBS); roots (CR); soil organic matter (CSOM); litter layers (CLITT); and the C export of woody biomass (CEXP), relative to the available incoming C from gross photosynthesis (GPP). We thus define the carbon sequestration efficiency (CSE) as the ratio of either modelled 5-year NEP or modelled lifetime NECB to modelled GPP in a given environment, constrained by climate, nitrogen availability and other factors included in the BASFOR model:

with

The modelled CSE_5-year can be contrasted with observation-based CSE_obs ($={\mathrm{NEP}}_{\mathrm{obs}}/{\mathrm{GPP}}_{\mathrm{obs}}$) derived from flux tower data over a similar, relatively short time period compared with a forest rotation (see Flechard et al., 2020). By extension, the CSE_lifetime indicator quantifies the efficiency of C sequestration processes by a managed forest system, reflecting not only biological and ecophysiological mechanisms, but also the long-term impact of human management through thinning frequency and severity.

For the N budget we define, by analogy to CSE, the N uptake efficiency (NUPE) as the ratio of N immobilized in the forest system to the available mineral N, i.e. the ratio of tree N uptake (N_upt) to the total N_supply from internal SOM mineralization and N-cycling processes (N_miner) and from external sources such as atmospheric N deposition (N_dep):

with

The fraction of N_supply not taken up in biomass and lost to the environment (N_loss) comprises dissolved inorganic N leaching (N_leach) and gaseous NO and N₂O emissions (N_emission):

Note that (i) NUPE is a different concept from the nitrogen use efficiency (NUE), often defined as the amount of biomass produced per unit of N taken up from the soil, or the ratio NPP∕N_upt (e.g. Finzi et al., 2007), and (ii) biological N₂ fixation, as well as N loss by total denitrification, are not accounted for in the current BASFOR version; also, leaching of dissolved organic N and C (DON, DOC) and dissolved inorganic C (DIC) is not included either, all of which potentially impact budget calculations.

2.2.5 Meta-modelling as a tool to standardize EC-based productivity data

One purpose of BASFOR modelling in this study was to gain knowledge on patterns of C and N fluxes, pools and internal cycling that were not, or could not be, evaluated solely on the basis of the available measurements (for example, SOM mineralization and soil N transfer; retranslocation processes at the canopy level; patterns over the lifetime of a stand). The model results were used to complement the flux tower observations to better constrain elemental budgets and assess the potential and limitations of C sequestration at the European forest sites considered here. Additionally, we used meta-modelling as an alternative to multivariate statistics (e.g. stepwise multiple regression, mixed non-linear models, residual analysis) to isolate the importance of N_r deposition from other drivers of productivity. This follows from the observations by Flechard et al. (2020) that (i) N_r deposition and climate were not independent in the dataset and that (ii) due to the large diversity of sites the limited size of the dataset and incomplete information on other important drivers (e.g. stand age, soil type, management), regression analyses were unable to untangle these climatic and other inter-relationships from the influence of N_r deposition.

BASFOR (or any other mechanistic model) is useful in this context, not so much to predict absolute fluxes and stocks but to investigate the relative importance of drivers, which is done by assessing changes in simulated quantities when model inputs are modified. Meta-modelling involves building and using surrogate models that can approximate results from more complicated simulation models; in this case we derived simplified relationships linking forest productivity to the impact of major drivers, which were then used to harmonize observations from different sites. For example, running BASFOR for a given site using meteorological input data from another site, or indeed from all other sites of the network, provides insight into the impact of climate on GPP or NEP, all other factors (soil, vegetation structure and age, N_r deposition) being equal. Within the boundaries of the network of 22 selected sites, this sensitivity analysis provides relative information as to which of the 22 meteorological datasets is most, or least, favourable to growth for this particular site. This can be repeated for all sites (22×22 climate scenario simulations). It can also be done for soil physical properties that affect the soil water holding capacity (texture, porosity, rooting depth), in which case the result is a relative ranking within the network of the different soils for their capacity to sustain an adequate water supply for tree growth. The procedure for the normalization of data between sites is described hereafter.

Additional nitrogen affects C uptake primarily through releasing N limitations at the leaf level for photosynthesis (Wortman et al., 2012; Fleischer et al., 2013), which scales up to GPP at the ecosystem level. Other major factors affecting carbon uptake are related to climate (photosynthetically active radiation, temperature, precipitation), soil (for example water holding capacity) or growth stage (tree age). In the following section, we postulate that observation-based gross primary productivity (GPP_obs), which represents an actuation of all limitations in the real world, can be transformed through meta-modelling into a standardized potential value (GPP^*) for a given set of environmental conditions (climate, soil, age) common to all sites, thereby enabling comparisons between sites. We define GPP^* as GPP_obs being modulated by one or several dimensionless factors (f_X):

where the standardization factors f_CLIM, f_SOIL and f_AGE are derived from BASFOR model simulations corresponding to the CEIP–NEU time interval around 2005–2010, as described below. The factors involved in Eq. (10) address commonly considered drivers but not nitrogen, which is later assessed on the basis of GPP^* rather than GPP_obs. Other potentially important limitations such as non-N nutrients, soil fertility, air pollution (O₃), poor ecosystem health and soil acidification are not treated in BASFOR and cannot be quantified here. Further, the broad patterns of the GPP vs. N_dep relationships reported in Flechard et al. (2020), i.e. a non-linear increase and eventual saturation of GPP as N_dep increases beyond a critical threshold, did not show any marked difference between the three forest PFTs (deciduous, pine, spruce), possibly because the datasets were not large enough and fairly heterogeneous. Thus, although PFT may be expected to influence C–N interactions, we did not seek to standardize GPP with an additional f_PFT factor.

To determine the f_CLIM and f_SOIL factors, the model was run multiple times with all climate and soil scenarios for the n (=22) sites, a scenario being defined as using model input data or parameters from another site. Specifically, for f_CLIM, the model weather inputs at each site were substituted in turn by the climate data (daily air temperature, global radiation, rainfall, wind speed and relative humidity) from all other sites; and for f_SOIL, the field capacity and wilting point parameters (Φ_FC, Φ_WP) and soil depth that determine the soil water holding capacity at each site (SWHC = (Φ_FC−Φ_WP) × soil depth) were substituted in turn by parameters from all other sites. Values of f_CLIM and f_SOIL were calculated for each site in several steps, starting with the calculation of the ratios of modelled GPP from the scenarios to the baseline value GPP_base such that

where i(1…n) denotes the site being modelled and j(1…n) denotes the climate dataset (j_CLIM) or soil parameter set (j_SOIL) used in the scenario being simulated (see Table S2 for the calculation matrices). The value of the X(i,j) ratio indicates whether the jth scenario is more (>1) or less (<1) favourable to GPP for the ith forest site.

For each site, the aim of the f_CLIM factor (and similar reasoning for f_SOIL) (Eq. 10) is to quantify the extent to which GPP differs from a standard GPP^* value that would occur if all sites were placed under the same climatic conditions. Rather than choose the climate of one particular site to normalize to, which could bias the analysis, we normalize GPP to the equivalent of a mean climate, by averaging BASFOR results over all (22) climate scenarios (Eqs. 14–15). However, since each of the scenarios has a different mean impact across all sites ($\stackrel{\mathrm{‾}}{X\left(j\right)}$, Eq. 12), we first normalize X(i,j) to the $\stackrel{\mathrm{‾}}{X\left(j\right)}$ value within each jth scenario (Eq. 13):

The normalization of X(i,j) to X_norm(i,j) ensures that the relative impacts of each scenario on all n sites can be compared between scenarios. The final step is the averaging for each site of X_norm(i,j) values from all scenarios (either j_CLIM or j_SOIL) into the overall f_CLIM or f_SOIL values:

or

The factors f_AGE were determined by first normalizing modelled GPP (base run) to the value predicted at age 80 for every year of the simulated GPP time series at those m (=12) mature sites where stand age exceeded 80. The age of 80 was chosen since this was the mean stand age of the whole network. The following ratios were thus calculated:

where k(1…m) denotes the mature forest site being modelled. A mean temporal curve for f_AGE (normalized to 80 years) was calculated to be used subsequently for all sites using the following:

4.3.1 From nutrient limitation to nitrogen saturation

Can nutrient limitation (nitrogen or otherwise) impact ecosystem carbon sequestration efficiency? Soil fertility has been suggested to be a strong driver at least of the forest biomass production efficiency (BPE), defined by Vicca et al. (2012) as the ratio of biomass production to GPP, with BPE increasing in their global dataset of 49 forests from 42 % to 58 % in soils with low to high nutrient availability, respectively. The study by Fernández-Martínez et al. (2014) of 92 forest sites around the globe reported a large variability in CSE ($=\mathrm{NEP}/\mathrm{GPP}$ calculated from FLUXNET flux data), which they suggest is strongly driven by ecosystem nutrient availability (ENA), with CSE levels below 10 % in nutrient-poor forests and above 30 % in nutrient-rich forests. The range of CSE values derived from flux measurements in our study (CSE_obs in Table 2) was similarly large, even though all of our sites were European and our dataset size was one-third of theirs (N=31, with 26 sites in common with Fernández-Martínez et al., 2014). We did not attempt in this study to characterize a general indicator of ENA beyond total N_r deposition; but if we use the high, medium or low (H, M, L) scores of ENA attributed to each site through factor analysis of nutrient indicators by Fernández-Martínez et al. (2014), we find that the H group (7 sites) has a mean CSE_obs of 32 % (range 16 %–48 %), the M group is slightly higher (7 sites, mean 39 %, range 21 %–61 %) and the L group has indeed a significantly smaller mean CSE_obs of 14 % (12 sites, range −9 % to 38 %). Interestingly, the mean N_dep levels for each group are H = 1.5 g N m⁻² yr⁻¹ (range 0.5–2.3 g N m⁻² yr⁻¹), M = 2.1 g N m⁻² yr⁻¹ (range 1.1–4.2 g N m⁻² yr⁻¹) and L = 1.3 g N m⁻² yr⁻¹ (range 0.3–4.1 g N m⁻² yr⁻¹); i.e. the highest mean CSE_obs of the three groups is found in the group with the highest mean N_dep (M).

The nutrients and other indicators of fertility considered by Fernández-Martínez et al. (2014) included, in addition to N, P, soil pH, C∕N ratios and cation exchange capacity, as well as soil texture and soil type. However, very few sites were fully documented (see their Supplement Table S1), data were often qualitative and other key nutrients were not included in the analysis (K, Mg and other cations; S also has been suggested to have become a limiting factor for forest growth following emission reductions; see Fernández-Martínez et al., 2017). The extent to which the overall fertility indicator quantified by ENA was driven by nitrogen in the Fernández-Martínez et al. (2014) factor analysis is not evident. At sites where other nutrients are limiting, the response to N additions would be small or negligible regardless of whether N itself is limiting. This places severe constraints on the interpretation of productivity data in response to N_dep, since most current models, which do not account for other nutrient limitations, cannot be called upon to normalize for differences between sites.

The impact of the fertility classification on CSE of the sites included in Fernández-Martínez et al. (2014) was questioned by Kutsch and Kolari (2015) on the basis of unequal quality of the EC flux datasets found in FLUXNET and other databases. By excluding complex terrain sites (and young forests) from the Fernández-Martínez et al. (2014) dataset, Kutsch and Kolari (2015) calculated a much reduced variability in CSE, with a reasonable mean value of 15 % (range 0 %–30 %), suggesting a much lower influence of nutrient status than claimed by Fernández-Martínez et al. (2014). In their reply, Fernández-Martínez et al. (2015) reanalysed the same subset of sites selected by Kutsch and Kolari (2015) but using the same generalized linear model as used in their original analysis of the whole dataset as opposed to the linear model used by Kutsch and Kolari (2015). Fernández-Martínez et al. (2015) then maintained that the findings of the original study were still valid for the restricted dataset, i.e. that the nutrient status had a significant influence on CSE.

The smaller European dataset of our study poses a similar dilemma. The much wider variation in CSE_obs than modelled CSE_5-year may both point to possible measurement issues if CSE_obs values (especially the larger ones) are considered ecologically implausible and/or inform on important ecological processes that are not accounted for in the model. Among the forests in our study that seemed particularly inefficient (CSE_obs<10 %) at retaining photosynthesized carbon (EN4, EN6, EN8, EN11, EN17, EB5), all were classified as L (low ENA) in Fernández-Martínez et al. (2014) and two (EN6, EN11) were even net C sources (R_eco>GPP). The EN4, EN6, EN17 sites had the three largest soil organic contents (SOCs, Fig. 10a), which may either have induced larger rates of heterotrophic respiration or may instead indicate low-fertility wet soils where both assimilation and respiration are suppressed. However, EN4 has also been reported as having unrealistically large ecosystem respiration rates (Anthoni et al., 2004). The EN8 site (mature pine-dominated forest in Belgium) was very unlikely to be N- or S-limited, having been under the high-deposition footprint of the Antwerp petrochemical harbour and local intensive agriculture for decades, even if emissions have declined over the last 20 years (Neirynck et al., 2007, 2011). However, the comparatively low LAI, GPP and CSE (Fig. 4 in Flechard et al., 2020) at this site are likely not independent of the historical N- and S-induced soil acidification, which has worsened the already low P and Mg availabilities (Janssens et al., 1999) and from which the forest is only slowly recovering (Neirynck et al., 2002; Holmberg et al., 2018). This site is actually an excellent example to illustrate the complex web of biogeochemical and ecological interactions, which further complicate the quantification of the (single-factor) N_dep impact on C fluxes. By not accounting for the low Mg and P availabilities and the poor ecosystem health, the BASFOR model massively overestimated GPP, R_eco and NEP at EN8 (Fig. 6 in Flechard et al., 2020). In fact, based on prior knowledge of this site's acidification history, and since such mechanisms and impacts are not mathematically represented in BASFOR, EN8 was from the start discarded from the calibration dataset for the Bayesian procedure (Cameron et al., 2018). The four lowest CSE_obs values were found at sites with topsoil pH < 4 (Fig. 10c), although other forests growing on acidic soils had reasonably large CSE_obs ratios.

The large variability in CSE_obs cannot be explained by any single edaphic factor (Fig. 10a–c), more likely by a combination of many factors that may include N_dep (Fig. 10e). As noted previously, C flux measurements at all four forest sites with N_dep>2.5 g N m⁻² yr⁻¹ (EN2, EN8, EN15, EN16) indicated lower productivity estimates than those in the intermediate N_dep range, or at least smaller than might have been expected from a linear N fertilization effect (Fig. 4 in Flechard et al., 2020). EN2 (spruce forest in southern Germany) is also well-documented as an N-saturated spruce forest with large total N losses (∼3 g N m⁻² yr⁻¹) as NO, N₂O and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Kreutzer et al., 2009), but its productivity and CSE are not affected to the same extent as EN8. Not all the difference is necessarily attributable to the deleterious impacts of excess N_r deposition, as suggested by the GPP normalization exercise (Fig. 8). For example, EN15 and EN16, planted on sandy soils, appear from meta-modelling to suffer from water stress comparatively more than the average of all sites (Fig. 6-Soil), if indicators of soil water retention based on estimates of soil depth, field capacity and wilting point can be considered reliable.

4.3.2 Forest age

Forest age is expected to affect photosynthesis (GPP), growth (NPP), carbon sequestration (NEP) and CSE for many reasons. A traditional view of the effect of stand age on forest NPP (Odum, 1969) postulated that R_aut increases with age and eventually nearly balances a stabilized GPP, such that NPP approaches zero upon reaching a dynamic steady state. Revisiting the paradigm, Tang et al. (2014) found that NPP did decrease with age (>100 years) in boreal and temperate forests, but the reason was that both GPP and R_aut declined, with the reduction in forest growth being primarily driven by GPP, which decreased more rapidly with age than R_aut after 100 years. However, the ratio NPP∕GPP remained approximately constant within each biome.

The effect of age on NEP and CSE is even more complex since this involves not only changing successional patterns of GPP and R_aut, but also those of R_het over a stand rotation of typically one century or more, which is much longer than the longest available flux datasets. Therefore age effects are often studied by comparing differently aged forest sites across the world, which introduces many additional factors of variation, including differences in water availability; soil fertility; or even tree species, genera, or PFTs. Forest and tree ages should in theory be normalized to account for species-specific ontogeny patterns; i.e. the age of 80 years may be relatively young for some species and quite old for others, and therefore population dynamics may be very different for the same age. Nevertheless, forest age has been suggested to be a dominant factor controlling the spatial and temporal variability in forest NEP at the global scale, compared with abiotic factors such as climate, soil characteristics and nutrient availability (Besnard et al., 2018). In that study, the multivariate statistical model of NEP, using data from 126 forest eddy-covariance flux sites worldwide, postulated a non-linear empirical relationship of NEP to age, adapted from Amiro et al. (2010), whereby NEP was negative (a net C source) for only a few years after forest establishment and then increased sharply above 0 (a net C sink), stabilized after around 30 years and remained at that level thereafter for mature forests (>100 years). This model, therefore, did not assume any significant reduction in forest net productivity after maturity, up to 300 years, consistent with several synthesis studies that have reported significant NEP of centuries-old forest stands (Buchmann and Schulze, 1999; Kolari et al., 2004; Luyssaert et al., 2008).

By analogy, our approach for accounting for the age effect was based on the modelled time course of GPP (Eqs. 16–17), which in the BASFOR model tended to stabilize after 100 years, and subsequently used a mean CSE that did not depend on stand age. However, the variability in CSE_obs appeared to be much larger in mature forests (>80 years) than in the younger stands (Fig. 10d). For the younger forests (<60 years, all sites probably still in an aggrading phase), the CSE_obs values were in a narrow band of 15 %–30 % and were well represented by model simulations, with the exceptions of EN1 and EB3 at around 50 % and of EN4 being near 0 % (with all three locations being high-elevation sites with complex terrain and potential EC measurement issues; see Flechard et al., 2020). By contrast, values for mature forests were either below 15 % or above 30 %. For some cold sites such as EN6 and EN11, growing in low-nutrient environments (e.g. peat at EN6) with high SOC (Fig. 10a) and/or high soil C∕N ratio (Fig. 10b) and low soil pH (Fig. 10c), or for the N-saturated and acidified EN8 site, the low CSE is not necessarily linked to age. Ageing, senescence and acidification may at some point curb sequestration efficiency in older forests, but even excluding the complex terrain sites, there remain a good number of productive mature sites with CSE_obs in the range 30 %–40 %, which questions the Odum (1969) paradigm of declining net productivity and C equilibrium in old forests.

4.3.3 Does nitrogen deposition impact soil respiration?

The overall net effect of N_r deposition on carbon sequestration must include not only productivity gains, but also indirect, positive or negative impacts on soil C losses, which all affect CSE. Carbon sequestration efficiency reflects the combined magnitudes of soil heterotrophic (R_het) and autotrophic (R_aut, both below- and above-ground) respiration components, relative to GPP. We postulated that the primary effect of N_dep and N_supply is on GPP, but potential side effects of N_dep or N additions on ecosystem and soil carbon cycling have been postulated. The traditional theory of the role of N on microbial decomposition of SOM was that, above a certain C∕N threshold value, the lack of N inhibits microbial activity compared with lower C∕N ratios (Alexander, 1977). However, reviews by Fog (1988) and Berg and Matzner (1997) found that microbial activity was often unaffected, or even negatively affected, by the addition of N to low-N decomposing organic material. The negative effects were mostly found for recalcitrant organic matter (high lignin content) with a high C∕N ratio (e.g. wood or straw), while N addition to easily degradable organic matter with a low C∕N ratio (e.g. leaf litter with low lignin content) actually boosted microbial activity. The meta-analysis by Janssens et al. (2010) of N manipulation experiments in forests suggests that excess N_r deposition reduces soil – especially heterotrophic – respiration in many temperate forests. They argue that the mechanisms include (i) a decrease in below-ground C allocation and the resulting root respiration, permitted by a lesser need to develop the rooting system when more N is available (see also Alberti et al., 2015); (ii) a reduction in the activity, diversity and biomass of rhizospheric mycorrhizal communities (see also Treseder, 2008); (iii) a reduction in the priming effect, the stimulation of SOM decomposition by saprotrophic organisms through root and mycorrhizal release of energy-rich organic compounds; (iv) N-induced shifts in saprotrophic microbial communities, leading to reduced saprotrophic respiration; and (v) increased chemical stabilization of SOM into more recalcitrant compounds. The authors point out that in N-saturated forests different processes and adverse effects are at play (e.g. base cation leaching and soil acidification). Of the five aforementioned mechanisms potentially involved in the suppression of soil respiration by N addition, only the first one (control by N availability of the root/shoot allocation ratio) is functional in BASFOR, and therefore our simulations do not include the other inhibitory effects of excess N on mycorrhizal, fungal and bacterial respiration.

An important implication of the negative impact of N_r on soil respiration is that the nitrogen fertilization effect on gross photosynthesis would be roughly doubled, in terms of C sequestration, by the concomitant decrease in soil respiration. In their meta-analyses of N addition experiments in forests and comparison of sites exposed to low vs. elevated N_dep, Janssens et al. (2010) show that both R_het and soil carbon efflux (SCE), a proxy for total R_soil ($=R_{\mathrm{het}}+R_{\mathrm{aut},\mathrm{soil}}$), tend to decline with N addition, be it through fertilization or atmospheric deposition, although the effect is far from universal. The negative N_dep response of R_het was much more pronounced for SOM than for leaf litter and stronger at highly productive sites than at less productive sites. The negative impact on SCE was mostly found at sites where N was not limiting for photosynthesis. When N is strongly limiting, and in young forests, N_r deposition may well favour SOM decomposition.

To examine the potential impact of N_dep on R_soil, we compiled the soil respiration data available from the literature and databases for the collection of forest sites in our study, which covers the whole N limitation to N saturation spectrum. Sites ranged from highly N limited boreal systems, where an N addition might trigger enhanced tree growth, increased microbial biomass and heterotrophic respiration, to N-saturated, acidified systems (EN2, EN8, possibly also EN15, EN16), in which poor ecosystem and soil health may lead to different ecological responses than those of the below-ground carbon-cycling scheme in Janssens et al. (2010).

Since the below-ground autotrophic (root and rhizosphere) respiration component is regulated to a large extent by photosynthetic activity (Collalti and Prentice, 2019), as well as seasonality in below-ground C allocation (Högberg et al., 2010), and contributes a large part of R_soil on an annual basis (Korhonen et al., 2009), the relationship of R_soil to N_dep is examined by first normalizing to GPP (Fig. 11a), yielding a soil respiration metric that is comparable between sites (for R_soil data, see Table S7 in the Supplement to Flechard et al., 2020). Similarly, the ratio R_soil∕R_eco shows the relative contribution of below-ground to total (ecosystem) respiration (Fig. 11c). Note that caution is needed when considering both R_soil∕GPP and R_soil∕R_eco ratios, since significant uncertainty may arise from (i) methodological flaws in comparing chamber versus eddy covariance measurements (e.g. considerations over tower footprint, spatial heterogeneity and representativeness of soil collars), (ii) uncertainty in deriving GPP and R_eco estimates from EC-NEE measurements, and (iii) different time spans for the EC and soil chamber measurements, affected by inter-annual flux variability. Thus, values of R_soil∕R_eco above unity (Fig. 11c), although physically nonsensical, do not necessarily imply large measurement errors but possibly also that there may be no spatial or temporal coherence in EC and chamber data (Luyssaert et al., 2009).

Figure 11Variability of normalized soil respiration metrics as a function of nitrogen deposition (a, c, e) and soil organic carbon (b, d, f). In all plots, the colour scale indicates mean annual temperature (MAT), and the symbol size is proportional to mean annual precipitation (MAP). R_soil: total soil respiration; R_eco: total ecosystem respiration; R_het: heterotrophic component of R_soil; GPP: gross primary productivity. DB: deciduous broadleaf; EN: coniferous evergreen needleleaf; EB: Mediterranean evergreen broadleaf; MF: mixed needleleaf–broadleaf forests.

Download

Either ignoring such outliers or judging that a measurement bias by soil chambers affects all sites the same way (e.g. systematic overestimation of soil respiration in low-turbulence conditions when using static chambers, Brændholt et al., 2017), we may argue that the apparent decrease in both chamber ∕ EC ratios R_soil ∕ GPP and R_soil ∕ R_eco with N_dep (Fig. 11a, c) has some reality, even if their absolute values are biased. Soil CO₂ efflux tends to be a larger fraction of GPP (>0.5) at the smaller N_dep rates (<1.5 g N m⁻² yr⁻¹) than at sites with larger N_dep, where this fraction is more often in the range 0.4–0.5. It is also noteworthy that the largest R_soil ∕ GPP ratios (EN5, EN17) are found at sites with relatively large SOC and topsoil C∕N ratio compared with the other sites (Figs. 10a and 11b). The R_het ∕ R_soil ratio also tends to decrease with N_dep (Fig. 11e), and although measured by different methods at the different sites, this is arguably a more robust metric than chamber ∕ EC respiration ratios, because the differential respiration measurements on control and treatment plots (root exclusion, trenching, girdling) are made on the same spatial and temporal scales.

Many other factors that impact soil respiration (age, soil pH, microbial abundance and diversity, etc.) are not considered here and are beyond the scope of this paper. In view of these uncertainties, if the assessment within this restricted dataset does not provide full and incontrovertible proof of the negative impact of N_r deposition on soil respiration, it at least is not in open contradiction to the prevailing paradigm that both below-ground autotrophic and heterotrophic respiration are expected to decrease as N_r deposition increases. However, the decreasing trends observed in Fig. 11a, c and e are largely driven by these few high-N_dep sites (>3 g N m⁻² yr⁻¹) in which the negative effects of N saturation and acidification very likely outweigh the benefits of reduced soil respiration in terms of C sequestration.