Nitrogen deposition in forests: Statistical modeling of total deposition from throughfall loads

1 Introduction

Excess atmospheric deposition of nitrogen, either in wet, dry or occult form, can result in a cascade of environmental effects impacting ecosystems diversity and functioning. While the wet deposition flux (WD) of nitrogen (N) is easily measured, the dry deposition flux (DD) is an important but very uncertain component in total N deposition flux (TD) estimations. The determination of the TD from the throughfall deposition flux measurements (TF) requires the quantification of canopy exchange (Draaijers and Erisman, 1995). There are large conceptually weaknesses and uncertainties in the application of canopy budget models for the canopy exchange quantification (Ahrends et al., 2020). Other methods applied are micrometeorological measurement and inferential modeling (IFM). The latter method is an approach by which deposition rates can be estimated on a routine basis, using measured concentrations of gasses and particles in the ambient air at a defined reference height and site- and element-specific deposition velocities (Hicks et al., 1991; Wesely and Hicks, 2000). A variety of approaches is used to determine deposition velocities (v_d), ranging from dynamic models based on stomatal conductance and atmospheric conditions to land-use specific empirical long-term averaged v_d (e.g., Schrader and Brümmer, 2014; Bytnerowicz et al., 2015; García-Gómez et al., 2018). At an intermediate level of complexity, published v_d are adapted for site specific conditions based on semi-empirical correction factors (Schmitt et al., 2005; Kirchner et al., 2014).

In addition, DD of N to forests may be estimated based on measurements with surrogate surfaces, such as e.g., Teflon strings placed under of roof, in combination with net throughfall estimates of an element that is regarded as biologically inert, such as e.g., sodium (Ferm and Hultberg, 1999; Karlsson et al., 2019). Following Lindberg and Owens (1992) the net throughfall flux (netTF) is the difference between throughfall deposition flux (TF) and WD (netTF = TF-WD) and represent the DD of an inert element to the forest canopies. Concentration ratios to sodium from the surrogate surface samples then give the DD of the element of interest in relation to the DD of sodium. This approach was tested for a wide range of N DD levels across Sweden (Karlsson et al., 2019, 2022).

For assessing the effect of N deposition on forest ecosystems, nitrogen gradient studies can provide an important knowledge contribution. These gradients are either based on modeled TD with varying spatial grid sizes (see the discussion on this issue in Braun et al., 2017) or on TF measured on the site. For example, TF and open field deposition measurements (together with stemflow in beech stands) are the mandatory methods for plots operated under the UNECE International Co-operative Program on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests) (ICP Forests, 2022). Measuring TF is a cheap, simple and low-maintenance method. Accordingly, the number of sites with long-term TF measurements is many times higher than the number of sites with air concentration measurements for applying the IF method (e.g., Ahrends et al., 2020, supplement). Therefore, many gradient studies rely on the results of TF. TF consists of the precipitation collected below canopy using open funnel-type samplers (McDowell et al., 2020). In parallel to throughfall measurements, precipitation is collected in the open field. Mostly, continuously open funnel-type samplers are used, and the open-field deposition is then called bulk deposition (BD). The BD samplers collect mainly WD with a small contribution of DD to the sampling equipment. In trees with a smooth bark like beech, stemflow (SF) is sampled, too, which is the substance flux with the water running down the stem. It contains similar components of inorganic ions as throughfall and can make up 5–20% of the sum of TF and SF (Neal et al., 1993; Chang and Matzner, 2000; Staelens et al., 2007) but it has not been included in gradient studies (e.g. van der Linde et al., 2018 which used the data of Waldner et al. (2014). Stemflow is less important for coniferous tree species, in particular trees with downward pointing branches, such as Norway spruce (Moldan et al., 1998; Kovács and Horváth, 2004; Małek, 2010). Similarly, in broadleaf evergreen trees, stemflow only contributes 3% to total rainfall (Rodrigo et al., 2003).

Throughfall collects WD, occult deposition (OD) and the part of the DD to the forest canopy, but canopy processes can also affect the deposition estimation. The fraction of DD collected varies depending on the atmospheric species and canopy exchange processes (Fenn et al., 2015). Part of the deposited N can be taken up by the forest canopy (Draaijers et al., 1996; Harrison et al., 2000; Adriaenssens et al., 2012; Karlsson et al., 2019), leading to an underestimation of TD by TF. To estimate the N fraction taken up by the canopy, canopy budget models can be applied which were developed by Ulrich (1994) and modified many times (Draaijers and Erisman, 1995; De Vries et al., 2001; Staelens et al., 2008). However, these models cannot estimate the direct uptake of gaseous compounds through the stomata (NH₃, NO₂) or the cuticle (NH₃). Moreover, canopy budget models might not be suitable to represent canopy processes in forests growing in semi-arid environments such as in the Mediterranean biogeographical region (Aguillaume et al., 2017). Further complications arise from transformation into organic N, which is also part of the canopy exchange process (Zimmermann et al., 2006; Cape et al., 2010; Massad et al., 2010; Adriaenssens et al., 2012) or from transformation of reduced forms of N into nitrate by epiphytic nitrifiers (Guerrieri et al., 2020). Other canopy processes such as NH₃ emission from leaves (Massad et al., 2010) or atmospheric N fixation by phyllosphere microorganisms (Rico et al., 2014) represent further inputs of N into the canopy whose possible contribution to TF N enrichment needs further research.

At the end, DD calculated with the different approaches needs to be combined with estimations of WD and OD to yield TD. Including OD, i.e., the deposition via cloud or fog water, is particularly relevant in mountain forests which are merged in clouds during a substantial part of the year (Bridges et al., 2002; Blaś et al., 2012; Hůnová et al., 2021a,b).

TF based N gradient studies give valuable information on N effects in ecosystems (e.g., van der Linde et al., 2018). To set critical loads of N, however, information on TD is needed. In the last years various studies were conducted to compare the measurement of N in TF to estimations of TD (García-Gómez et al., 2018; Karlsson et al., 2019; Thimonier et al., 2019; Ahrends et al., 2020). This allows finding relations between TF and TD and deriving a transfer function. The aim of this study is to provide TD estimates for forests where only TF of N is available.

2 Materials and methods

Data were collected from studies which measured total N deposition either using the IFM for DD estimates (Spain: García-Gómez et al., 2018; Switzerland: Thimonier et al., 2019; Germany: Ahrends et al., 2020) or an equivalent estimate using teflon strings (Sweden: Karlsson et al. (2019), along with BD (or WD) and TF. The Thimonier et al. (2019) dataset was amended with TF from an additional site from Braun et al. (2018). Additional datasets included were from Germany (Zimmermann et al., 2006) and Belgium (Neirynck et al., 2007), various NitroEurope sites across Europe (Flechard et al., 2011), Finland (Korhonen et al., 2013) as well as 36 plots from the Czech Republic (Hůnová et al., 2016). The final dataset consisted of 103 coniferous, 33 deciduous, 4 broadleaved evergreen and 2 mixed stands with data from 1–15 years each. The distribution of the plots is shown in Figure 1, a description of the datasets is given in Supplementary Table 1. Since stemflow has not been considered in many gradient studies (Waldner et al., 2014), it was not included even when it was measured.

Figure 1. Location of the plots included in the current study. Broadleaf deciduous, broadleaf evergreen, and mixed were combined to one forest group in the analysis.

Response variable of the regression analysis was total N deposition. Explanatory variables were included according to theoretical considerations and availability. The initial model included covariates such as precipitation, altitude, annual temperature, longitude, forest type, year of the measurement and, to cover for differences in measurement methodology, also country. We calculated a Variance Inflating Factor (VIF) to check for collinearity. The data were subjected to an explorative analysis using mixed regression, function glmer in R, and polynomial functions for non-linearity (Bates et al., 2015). Predictors increasing the Akaike information criterion, AIC (Burnham and Anderson, 2002) were removed, and transformations of TD or TF were tested to yield a uniform and normal distribution of the residuals. The model resulting from this selection procedure used log_e transformation for both TD and TF, forest type as well as year of measurement and country as random variables. Forest type was included as binary variable (0 = coniferous forests, 1 = broadleaf deciduous, broadleaf evergreen and mixed forests). The final multivariate mixed model (MMM) has the form:

However, during model setup, it was noticed that the inclusion of covariates led to mismatch of the intercept with the data, resulting in systematically overestimation of low values. The final selected model had therefore the form:

without any random effects and covariates. The residuals were checked for outliers, homoscedasticity and normal distribution. However, the model with covariates is listed for information, too.

3 Results

Table 1 lists the output from the MMM model. Forest type (coniferous vs. broadleaved deciduous, broadleaved evergreen and mixed) and geographical longitude were significant predictors. The VIF analysis suggested, however, that longitude was too strongly correlated with TF so it was removed. The coefficient for forest type suggested that the estimates for TD are lower for deciduous stands than for coniferous stands. The model in Table 1 resulted in an estimate of the intercept which was higher than the majority of the data points (Figure 2). Therefore, the model was simplified to the one shown in Table 2, which considerably improved the intercept estimates (Figure 3). Because of the many points from Sweden with low TD, the intercept could be calculated rather precisely. As an example, the estimate for 2 kg N ha^–1 yr^–1 TF is 6 kg N ha^–1 yr^–1 TD. For 5 kg N ha^–1 yr^–1 TF, the estimate for TD is on average 10 kg N ha^–1 yr^–1.

Table 1. Regression output of the multivariate mixed model (MMM).

Figure 2. Output of the multivariate mixed model: raw data plot and regression with 95%-confidence intervals as gray fields. Points are semi-transparent, so darker dots indicate overlap of several points.

Table 2. Regression output of the final selected model.

Figure 3. Output of the univariate model: raw data plot and regression with 95%-confidence intervals as gray fields. Points are semi-transparent, so darker dots indicate overlap of several points.

4 Discussion

The measurement of TD is complicated as there is a number of components contributing to it with varying importance for different ecosystems. The IFM for measuring TD uses a number of assumptions on v_d but is well established as reference method. In some plots it has been validated with more complicated eddy flux measurements (Neirynck et al., 2007). The teflon string method used in Scandinavia relies on another principle but seems to yield comparable results (Karlsson et al., 2019).

The results of this study clearly confirm that canopy uptake occurs, and that it is significant. The resulting shape, with values above the 1:1 line for TF < 20 kg N ha^–1 yr^–1 (Figure 4), might be explained by gradients of N deposition and demand in Europe. In far northern latitudes, both the demand and the supply are low, hence the difference between TF and TD will be small in absolute terms (Karlsson et al., 2019). At more central latitudes in Europe, the conditions for physiological activity are better, while the N deposition covers a large range from deposition hotspots to less polluted areas. At some of the less polluted sites, N demand is high while N is not available in sufficient quantities, resulting in a large canopy uptake rate and, correspondingly, in a large difference between TF and TD (with TF < TD). The difference between TF and TD is largest between 4 and 8 kg N ha^–1 yr^–1 in TF (see also Karlsson et al., 2019). At these more polluted sites the uptake will be low, and hence the difference will be low. In Mediterranean areas, as published in Meixner and Fenn (2004) and García-Gomez et al. (2016), canopy uptake and retention strongly varies throughout the year. During dry periods (mainly summer, but also winter), when climatic conditions determine low physiological activity and N uptake, there is an accumulation of deposited atmospheric N, resulting in N-enriched TF. On the other hand, during spring, there is a higher N demand, reducing the content of N in TF. This asynchrony between N deposition and demand in Mediterranean ecosystems could cause a strong variability on the relationship between TF and TD. Since this regression has been obtained from a large geographical range of studies within Europe, it is meant to be a generalized tool to estimate TD in cases where only TF is measured.

Figure 4. Total N deposition regressed against throughfall N deposition after back-transformation to linear scales (univariate model). The solid line shows the regression. gray range 95% confidence interval. Dashed line: 1:1 line. Blue field: Range where the regression should not be used (TF > TD). BE, Belgium; DK, Denmark; NL, Netherlands; CH, Switzerland; CZ, Czech Republic; DE, Germany; ES, Spain; FI, Finland; SE, Sweden.

The fitted relation of TD vs. TF indicates that TF is lower than TD at TF rates above around 20 kg N ha^–1 yr^–1. To our knowledge, the leaching of inorganic nitrogen from forest canopies has not been reported in the literature until now. On the contrary, numerous reports indicate that inorganic nitrogen may be taken up by canopy foliage, epiphytic lichens or other microflora (Sparks, 2009; Butterbach-Bahl et al., 2011; Woods et al., 2012; Guerrieri et al., 2021). In view of the large scatter of modeled TD at high TF rates, we hypothesize that the data points with TD < TF (i.e., data points below the 1:1 line in Figure 4, at TF values > 20 kg N ha^–1 yr^–1) reflect the uncertainty inherent in both the TF measurements and the TD estimations. It is thus not advisable to use the relation reported in our study at TF rates beyond 20 kg N ha^–1 yr^–1. Instead, TF should be used as a lower boundary estimate of TD in regions where our relation indicates TD < TF. This recommendation is motivated by the realization that the analyzed dataset contains a wide cascade of uncertainties that can accumulate with each measurement and/or calculation step of deposition estimation:

(A) Uncertainty in throughfall sampling and missing stemflow: The stand deposition flux (SD) is the sum of throughfall and stemflow and would be the better basis for comparison with TD. However, stemflow is not always provided in the ICP Forests measurements, and gradient studies usually work with TF only, thus introducing a systematic underestimation of stand deposition in deciduous stands with a smooth bark (e.g., beech). However, this mainly concerns only the beech tree species (Clarke et al., 2016). As in the present study the data analysis had to be restricted to throughfall, the term “throughfall” is used in this study even when “stand deposition” would be more appropriate. At dry sites, precipitation might be insufficient to mobilize the N compounds intercepted by the canopy during dry periods, which leads to lower throughfall fluxes than theoretically expected (Bytnerowicz et al., 2015; Thimonier et al., 2019).

The ICP Forests manual (Clarke et al., 2016) suggests a number of samplers depending on forest stand heterogeneity such that the resulting uncertainty in throughfall measurements is 20% (Draaijers et al., 1998; Thimonier, 1998; Maes et al., 2014). For various reasons, such as occasional contamination by biological material in the samplers, etc., this might be rather a lower boundary estimate.

(B) Uncertainty in TD of those N deposition pathways covered by the methods. Flechard et al. (2011) demonstrated that discrepancies between N deposition estimates from different IFM at one site can be several times larger than between sites. One possible reason for this is that the IFM usually is a combination of observations and modeling with numerous assumptions. The gaseous components NO₂ and NH₃ were assessed at the intensive monitoring sites. The accuracy level of these measurements is stated to be about ± 30% (Schaub et al., 2016). Other aspects, such as the low sampling frequency of passive samplers, may cause biased results (Schrader et al., 2018). A probably larger uncertainty is harbored in the parameterization of v_d. Several mechanistic models for v_d exist, which e.g., utilize detailed information on meteorology and soil water availability to model stomatal opening. This allows for the parameterization of stomatal exchange if ambient air concentrations compared to foliar concentrations are known (Wichink Kruit et al., 2012). While the complexity of these comprehensive models offers the opportunity for an accurate quantification of DD fluxes, additional uncertainty arises from assumptions about input parameters. The studies used in this paper choose a different approach, starting from typical v_d from the literature. These forest type specific values were adjusted using the most important determinants for v_d in order to derive “semi-empirical correction factors.” However the different authors have used various v_d for the IFM: some such as Ahrends et al. (2020) used site-specific modifications on the v_d to account for differences in exposition, wind speed, etc., others differentiated the v_d with forest type (Thimonier et al., 2019) or simply used the annual mean v_d (Hůnová et al., 2016). Thus, the resulting v_d for gaseous and particulate substances are not easily comparable and can only be considered as a very rough estimate. Furthermore, the DD of HNO₃, particulate NH₄⁺ and particulate NO₃^– was often calculated from open field deposition based on an empirical relation, because local concentration measurements were rarely available. At the same time we are fully aware of the fact that v_d belongs to key factors in DD estimate, is highly variable due to its dependence on surface, meteorological conditions, season and time of the day (Brook et al., 1997; Zhang et al., 2009), and thus can introduce a substantial uncertainty in the deposition estimate.

(C) Some N deposition pathways not covered by the methods: Inorganic N deposition can be converted to organic N in the canopy by microorganisms (e.g., Cape et al., 2010). This N deposition pathway is accounted for in the TD estimates but is missing in the measured stand deposition of inorganic N. Similarly, the stomatal canopy uptake of NH₃ and NO₂ is accounted for in the TD estimates but not in the TF. The magnitude of these effects probably strongly varies depending on site-specific aspects such as plant species, meteorological conditions and soil water availability.

The final regression model gave a reasonable fit to the data and can serve to translate TF measurements into TD estimates in the desired context of gradient studies of forest ecosystem responses to N deposition. A number of covariates were tested in the regression model. Precipitation was not significant but forest type was a significant predictor in the mixed model, and the inclusion of country and year as random variables improved the model, too. Longitude was too closely related with throughfall deposition and confounded forest type. This is, however, without practical consequence for the final selected model as the inclusion of covariates affected the intercept estimates, resulting in a systematic overestimation of the N deposition at low throughfall values. The knowledge on a longitude effect would also not help to improve the result of a gradient study where detailed site information is not available anyway.

In their comparison of TD estimations from the EMEP chemical transport model (Simpson et al., 2020) with TF from ICP forest plots, Marchetto et al. (2021) came to a similar conclusion as the present study. They report that EMEP-modeled TD of NH₄⁺ is on average 77% higher and modeled TD of NO₃^– is on average 72% higher compared to TF measurements at the same site. Their scatter plots show the largest differences at NH₄⁺ loads in TF between 5 and 10 kg N ha^–1 yr^–1 and at NO₃^– loads in TF between 3 and 8 kg N ha^–1 yr^–1. As the EMEP estimate of total N deposition is now available at 0.1° grid size (which corresponds to between 3.5 km and 8.5 km in Europe), it can provide additional information for gradient studies. For example it can help to set limits to estimates from throughfall studies in low deposition areas. However, it has to be taken into account that in the EMEP grid various receptors in the grid cell are weighted. For the use with forests, the data can be downloaded with specific v_d for coniferous and deciduous forests (EMEP, 2022).

5 Conclusion

TF clearly underestimates TD of nitrogen. This underestimation is largest in areas with low to medium deposition. With the present data set covering large parts of Europe it is possible to derive a critical load estimate from gradient studies using TF data (Bobbink et al., 2022).

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number (s) can be found below: doi: 10.5061/dryad.2z34tmpqc.

Author contributions

SA and SB made the concept of the study. SB, BA, AS, and PK collected the data. RA, HG-G, IH, GK, and AT contributed to the data. All authors contributed to the article and approved the submitted version.

Funding

This evaluation was partly based on data that was collected by partners of the official UNECE ICP Forests Network (http://icp-forests.net/page/participating-countries). Part of the data was Co-financed by the European Commission. Partial funding of data collection was provided by the European Union Under Council Regulation (EEC) 3528/86 on the Protection of Forests Against Atmospheric Pollution, the Regulation (EC) 2152/2003 concerning monitoring of forests and environmental interactions in the community (Forest Focus) and by the project LIFE 07 ENV/D/000218 Further Development and Implementation of an EU-Level Forest Monitoring System (FutMon). The German part of evaluation was part of the project “Site adapted management in privately owned small forests” funded by BMEL/FNR (FKZ: 2220NR269C). The Swedish data was collected within the Swedish Throughfall Monitoring Network (SWETHRO). SWETHRO was run by IVL Swedish Environmental Research Institute and funded by Regional Air Quality Protection Associations, County Administrative Boards, and the Swedish Environmental Protection Agency. The Swiss study was partly funded by the Swiss Federal Office for the Environment (FOEN), the environmental departments of the cantons GR, AG, BE, GE, SO and TI, and OSTLUFT. The Spanish data were collected in the framework of the EDEN-Med Project (FEDER/MICINN–AEI CGL2017-84687-C2-1-R) and the collaboration in this study was co-funded by the Framework on atmospheric pollution and persistent organic pollutants between DGCEA and CIEMAT (ACTUA-MITERD). The Czech data used for the analysis were provided by the Czech Hydrometeorological Institute, Czech Geological Survey, and Forestry and Game Management Research Institute.

Acknowledgments

We thank Christian Schindler, Swiss Tropical and Public Health Institute, University of Basel, Switzerland, for statistical advice. We also thank Christophe Flechard, Johan Neirynck, Janne Korhonen, Klaus Butterbach-Bahl and Pasi Kolari for their assistance in locating additional data sets. We thank the numerous staff involved in operating the deposition monitoring in the different countries in the last decades. We acknowledge the Dry Deposition Group initiative within the Expert Panel on Deposition of ICP Forests (http://icp-forests.net/group/dry-deposition), which stimulated this collaboration.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2022.1062223/full#supplementary-material

References