Comparison of Methods for the Estimation of Total Inorganic Nitrogen Deposition to Forests in Germany

A reliable quantification of total inorganic nitrogen (TIN) deposition to forests is required for the evaluation of ecological effects of TIN inputs to forests and to monitor the success of clean-air policy. As direct measurements are scarce, different modeling approaches have been developed to estimate TIN deposition to forests. Three common methods are the (i) “canopy budget model,” (ii) “inferential method,” and (iii) “emission based estimates” using a chemical transport model. Previous studies have reported considerable and site-specific differences between these methods, complicating the interpretation of results. We use data from more than 100 German intensive forest monitoring sites over a period of 16 years for a cross-comparison of these approaches. Non-linear mixed-effect models were applied to evaluate how factors like meteorology, terrain and stand characteristics affect discrepancies between the model approaches. Taking into account the uncertainties in deposition estimates, there is a good agreement between the canopy budget and the inferential method when using semi-empirical correction factors for deposition velocity. Wet deposition estimates of the emission based approach were in good agreement with wet-only corrected bulk open field deposition measurements used by the other two approaches. High precipitation amounts partly explained remaining differences in wet deposition. Larger discrepancies were observed when dry deposition estimates are compared between the emissions based approach and the other two approaches, which appear to be related to a combination of meteorological conditions and tree species effects.

A reliable quantification of total inorganic nitrogen (TIN) deposition to forests is required for the evaluation of ecological effects of TIN inputs to forests and to monitor the success of clean-air policy. As direct measurements are scarce, different modeling approaches have been developed to estimate TIN deposition to forests. Three common methods are the (i) "canopy budget model," (ii) "inferential method," and (iii) "emission based estimates" using a chemical transport model. Previous studies have reported considerable and site-specific differences between these methods, complicating the interpretation of results. We use data from more than 100 German intensive forest monitoring sites over a period of 16 years for a cross-comparison of these approaches. Non-linear mixed-effect models were applied to evaluate how factors like meteorology, terrain and stand characteristics affect discrepancies between the model approaches. Taking into account the uncertainties in deposition estimates, there is a good agreement between the canopy budget and the inferential method when using semi-empirical correction factors for deposition velocity. Wet deposition estimates of the emission based approach were in good agreement with wet-only corrected bulk open field deposition measurements used by the other two approaches. High precipitation amounts partly explained remaining differences in wet deposition. Larger discrepancies were observed when dry deposition estimates are compared between the emissions based approach and the other two approaches, which appear to be related to a combination of meteorological conditions and tree species effects.

INTRODUCTION
During the last 70 years emissions of nitrogen (N) species to the atmosphere from traffic, industrial processes, and agriculture have drastically increased over pre-industrial levels and a significant decrease in the next decades in Europe is not expected (Simpson et al., 2014). The resulting atmospheric deposition of inorganic N to forests is an important determinant of tree growth (Etzold et al., 2020), rendering nitrogen deposition an essential input variable in decision support systems for forestry under environmental change (Panferov et al., 2011;Thiele et al., 2017) and carbon uptake (Du and De Vries, 2018). Accurate quantification of N deposition is also necessary for the estimation of nitrate leaching from forest ecosystems (MacDonald et al., 2002;Johnson et al., 2018;Vuorenmaa et al., 2018) and the calculation of N budget changes in forest soils (Fleck et al., 2019). On a political and administrative level, N deposition estimates are required to assess the success of clean air policy (Hettelingh et al., 2017), the exceedance of critical loads for eutrophication and acidification (De Vries et al., 2015) and in the context of licensing procedures for nitrogen emitting facilities.
The total deposition (TD) of N into forest ecosystems occurs via three pathways (Unsworth and Fowler, 1987): Wet deposition (WD) comprises deposition via rain, snow and hail; dry deposition (DD) consists of gases and particles deposited on surfaces or directly taken up by vegetation; and occult deposition (OD) refers to the deposition of fog. DD and OD to forests is typically larger compared to other land cover types, due to the large surface area of the canopy and their high aerodynamic roughness. The sum of DD and OD is also referred to as interception deposition (ID, Ulrich, 1994). Unlike WD, which is fairly easy to assess (Staelens et al., 2008;Dämmgen et al., 2013), the quantification of ID is much more challenging. As OD is often of orographic origin, it usually only contributes significantly to TD in mountainous regions (Kirchner et al., 2014;Hunová et al., 2016).
The accurate quantification of DD fluxes to forests is still challenging due to a large variety of N species, their chemical reactivity, a high uncertainty in the estimation of deposition velocities (Saylor et al., 2019) and different deposition pathways including bi-directional fluxes (Wichink Kruit et al., 2012). Currently, micrometeorological methods (e.g., eddy covariance and gradient techniques) are regarded as the most accurate approaches to quantify DD and OD (Marques et al., 2001;Mohr et al., 2005;Schmitt et al., 2005;Brümmer et al., 2020). However, micrometeorological methods require a considerable measurement effort and observational data are therefore typically only available for short observation periods at a limited number of locations.
Three other methods are frequently used where information on deposition of total inorganic nitrogen (TIN) is required. Firstly, the canopy budget model (CBM) approach developed by Ulrich (1994) and modified many times (e.g., Draaijers and Erisman, 1995;De Vries et al., 2001;Staelens et al., 2008) is applicable where assessments of open field and throughfall precipitation and element concentrations are available, e.g., for intensive forest monitoring sites (De Vries et al., 2003;Meesenburg et al., 2004Meesenburg et al., , 2016Talkner et al., 2010). The application of the Ulrich (1994) CBM is straightforward and requires no empirical parameters. However, estimates of TIN deposition with CBM are questionable due to debatable assumptions, limited understanding of canopy ion exchange processes and propagation of measurement errors in calculations (Staelens et al., 2008;Adriaenssens et al., 2013).
Secondly, for monitoring sites with observations of ambient air concentrations of major gaseous N species (e.g., NH 3 and NO 2 ), DD can be estimated using the inferential method (IFM), i.e., by multiplying the concentrations with deposition velocities (Zimmermann et al., 2006). A variety of approaches is used to inform deposition velocities, ranging from dynamic models based on stomatal conductance and atmospheric conditions to land-use specific empirical long-term averaged deposition velocities (e.g., Schrader and Brümmer, 2014). At an intermediate level of complexity, published deposition velocities are adapted for site specific conditions based on semi-empirical correction factors (Schmitt et al., 2005;Kirchner et al., 2014). The IFM approach suffers from considerable variances in deposition velocities as shown by different review studies (e.g., Staelens et al., 2012;Schrader and Brümmer, 2014). In addition, the observation of NO 2 and NH 3 ambient air concentrations yield an extra uncertainty of ±30% (Schaub et al., 2016). At the end, DD calculated with the IFM needs to be combined with measurements of WD to yield TIN TD.
The third method to estimate TIN deposition is a combination of emission inventories and a chemical transport model using meteorological data for the simulation of the regional circulation (referred to as emission based method, EBM, in the following). A range of modeling systems exist (Vivanco et al., 2018). For Germany, the German Environmental Agency has funded the development of an emission-based approach that yields features with a higher spatial resolution compared to some European scale models . Approaches with higher spatial resolution have been used in Germany for many years (Gauger et al., 2008;Builtjes et al., 2011;Schaap et al., 2015Schaap et al., , 2017Schaap et al., , 2018 and its results have been included in numerous impact studies (Hauck et al., 2012;Fleck et al., 2017Fleck et al., , 2019Thiele et al., 2017). In Germany, the EBM approach integrates emission inventories and a large number of local measurements of wet or bulk deposition to TIN deposition estimates with complete spatial coverage. Major challenges lie in the accuracy and spatiotemporal resolution of emission data and the parametrization of receptor-specific deposition processes with respect to DD and OD (Saylor et al., 2019).
Previous studies comparing TIN deposition derived with different methods occasionally reported considerable site-specific discrepancies. This is true for micrometeorological methods and CBM (Marques et al., 2001;Mohr et al., 2005), IFM and CBM (Schmitt et al., 2005;Kirchner et al., 2014) as well as EBM and CBM  or multiapproach intercomparisons (Brümmer et al., 2020). Contrarily, Zimmermann et al. (2006) found a good agreement between CBM and IFM, while in the study of Thimonier et al. (2019) a throughfall based method (although not a CBM), IFM and EBM generally yielded similar rates of TIN deposition with notable exceptions at some sites. These findings suggest that comparisons of deposition estimates are most informative when carried out across a large number of sites and long observation periods.
Model intercomparisons can support the interpretation of results from single methods and indicate conditions where Modeled based on assumptions about the concentration ratios of nitrogen compounds relative to Na + in WD, DD, and OD. Informed by measurements of bulk open field and stand precipitation In situ measurements of ambient air NH 3 and NO 2 concentrations using passive samplers combined with site specific deposition velocities. HNO 3 , NO − 3 , and NH + 4 estimated after Schmitt et al. (2005) Chemical transport model LOTOS-EUROS Taken from EBM Estimated based on modeled meteorological data, concentration fields used for the WD of EBM and empirical fog water enrichment factors specific approaches are more or less reliable. The aim of this study is to contribute a systematic comparison of three common approaches to estimate TIN deposition to forests for an extended geographic and temporal coverage. Based on data from around 100 intensive forest monitoring sites (42 sites for IFM) in Germany over a period of 16 years, we derive TIN deposition estimates with (i) the "canopy budget model" (CBM), (ii) the "inferential method" (IFM), and (iii) "emission based estimates" (EBM). We evaluate the discrepancies between the different approaches using mixed effect models in order to analyze if they are determined by spatial, temporal, meteorological, and site specific factors. Where applicable, model discrepancies were analyzed separately for wet deposition (WD), dry deposition (DD), and total deposition (TD). In this study, we hypothesize that the difference between the methods can be partly explained by meteorology, terrain characteristics, site specific factors and levels of ambient air concentrations.

MATERIALS AND METHODS
In the following subchapters a description of the assessments and models is provided. The overall study design is summarized in Table 1.

Study Sites and Data Coverage
In Germany, data collection under the International Cooperative Programme on Assessment and Monitoring of Air Pollutions Effects on Forests (ICP Forests) is conducted since more than two decades (Ferretti and Schaub, 2014). As part of this program, atmospheric deposition is assessed at intensive monitoring plots ("Level II") by means of precipitation sampling in the forest stands and at nearby open field sites (De Vries et al., 2003). Hundred and four sites with varying temporal coverage and a variety of forest stand types were examined for this study (Figure 1). (Larix decidua Mill.) and stands with more than one dominant tree species [13]. The plots are located at altitudes between 10 and 1,522 m a.s.l. Mean air temperature and mean annual precipitation (1981-2010) ranged from 2.6 to 10.9 • C and from 558 to 2,444 mm, respectively. While data according to EBM approach are available for all plots and years in the period 2000-2015, application of the CBM and the IFM approach is limited by the availability, completeness and quality of observations. The CBM and IFM were calculated for 1,237 and 194 site-years, respectively. Further information on the sites and data availability is provided in the supplementary material (Supplementary Table 1).

Sampling Procedures, Chemical Analysis and Data Quality
Deposition assessments for the CBM approach were conducted according to the ICP Forests Manual on sampling and analysis of deposition (Clarke et al., 2016). In short, open field deposition was collected by 3-6 continuously open bulk samplers at sites in the vicinity of the forest stands. Between 9 and 27 collectors are placed under the forest canopy in varying spatial arrangements in order to cover the spatial variation in throughfall deposition. At plots with European beech, stemflow is assessed at a subset of the trees (Clarke et al., 2016). Usually, samples from multiple samplers are pooled in order to reduce the analytical effort.
Samples are collected at least fortnightly, filtered, and then stored in the dark at about 4 • C or below before chemical analyses are performed. For some plots samples are mixed to monthly samples. Deposition samples were analyzed for Na + , NH + 4 , and NO − 3 , and a range of other parameters by different laboratories. Due to the standardized methods (Clarke et al., 2016), data are comparable and laboratory results are checked with currently recommended methods (Mosello et al., 2005): (1) the ion balance, (2) a comparison between measured and calculated conductivity, (3) a comparison between the sum of the inorganic forms of nitrogen and total nitrogen, and (4) the Na/Cl ratio. If analytical results are suspicious, analyses are repeated. The QA/QC procedures further included the use of control charts for internal reference material to check long-term FIGURE 1 | Location of 104 German intensive monitoring plots considered in the study. Symbols representing beech, oak, spruce and pine dominated forest sites and sites with other trees (Douglas fir, larch, mixed forest). comparability within national laboratories (König et al., 2013) as well as participation in periodic ring tests to check comparability between laboratories (Marchetto et al., 2006).
The three deposition fluxes (bulk open field, throughfall, stemflow) are then calculated by multiplication of the precipitation amount with the corresponding ion concentration in the analyzed precipitation samples and summed up to annual rates. Only annual open field and throughfall deposition fluxes, which were based on at least 292 days of collection (80% of the year) have been included. In case of data gaps (between 0 and 20% of the year), the corresponding deposition rate has been extrapolated to full temporal coverage based on the assumption that daily deposition fluxes in the unobserved time periods equal the average daily deposition fluxes in the observed period for the respective plot, sampler and year (Fischer et al., 2007). For stemflow (on average amounting to approximately 10-16% of the Na + , NH + 4 and NO − 3 throughfall deposition flux), longer data gaps frequently occurred. They were filled based on plotand substance-specific average ratios between stemflow and throughfall deposition rates for those measurement periods where throughfall was available. Finally, the stand precipitation (ST) deposition flux was calculated as the sum of throughfall and stemflow for beech plots and equaled the throughfall flux for all other plots.
Ambient air concentrations of NO 2 , and NH 3 were assessed at 43 sites by passive samplers in different years (Supplementary Table 1). Measurements were made using IVL passive samplers usually at 2 m above ground, installed in the open field where bulk open field precipitation was assessed (Schaub et al., 2016). The samples were mostly processed and analyzed by the Swedish Environment Research Institute (IVL). More details on the samplers are given in Swaans et al. (2007). All years with a data completeness of <80% for a specific air pollutant were excluded for the respective measurement site (Schaub et al., 2016).

Approaches to Estimate Total Inorganic Nitrogen Deposition
N deposition to forests consists of several N components. We refer to TIN as the substances NH + 4 and NO − 3 (ions in precipitation and aerosols) as well as NH 3 , NO 2 , and HNO 3 (gases). The three approaches compared in this study cover these components to a slightly different extent. The EBM reports deposition rates for oxidized N (NO y ) which also includes compounds like HNO 2 and NO. The canopy budget model also accounts for these substances as they mostly react to N forms captured by the measurements used for the method (Thimonier et al., 2019). The canopy budget model, however, only partly covers the deposition of gases in general (see below). The inferential method explicitly models the deposition fluxes of the five TIN compounds. These differences must be taken into account for the interpretation of results. The deposition of organic N is not the specific subject of this study. In the following, the three methods are briefly described.

Canopy Budget Model
A number of CBM versions exist and differences between models were evaluated in other studies (Staelens et al., 2008;Adriaenssens et al., 2013). We selected the approach of Ulrich (1994) as a relatively robust and conservative version [the approach probably underestimates TIN TD for several reasons (Meesenburg et al., 2009), see Discussion)]. Based on the assessment of NO − 3 and NH + 4 in bulk open field precipitation and ST, dry deposition of gaseous and particulate N species is estimated. Therefore, the sum of the calculated NO − 3 and NH + 4 TD is referred to as TIN TD. In detail, the CBM of Ulrich (1994) calculates the TD of the nitrogen components (NC) NO − 3 and NH + 4 as the sum of WD and interception deposition (ID): Wet deposition was estimated from bulk open field precipitation based on correction factors from Gauger et al. (2002). Due to the high temporal and spatial variability of the factors, this can only be considered as a rough approximation. More details on the limitations are given in section Uncertainties in Methods and Measurements. The ID is conceptually split into (1) particulate ID (ID part,NC ), consisting of particulate DD and OD, and (2) gaseous deposition (ID gas, NC ).
The model assumes that concentration ratios of substances (Na + :NC) in ID part are similar to concentration ratios of substances in WD. Furthermore, it is assumed that ID part of sodium can be estimated as the difference of ST and WD (zero net canopy exchange of Na + ). Based on these assumptions, ID part of each of the two nitrogen compounds can be estimated as: The ID of gaseous N species is estimated as the share of ST that is not explained by WD or ID part of NO − 3 and NH + 4 , respectively.
If the sum of WD and ID part exceeds ST, no gaseous deposition can be calculated (treated as zero). Finally TIN TD is calculated as the sum of TD of NH + 4 and NO − 3 . For more details, the underlying assumptions and limitations see Ulrich (1994) or Meesenburg et al. (2009). Note that the CBM yields an estimate of ID but does not allow to differentiate between DD and OD. Therefore, the OD and DD fluxes of the other two methods (IFM and EBM, see below) are also aggregated to yield the respective ID fluxes. As the OD share among the ID is usually relatively small, and in order to reduce the complexity of the figures and tables, we present the results for OD and DD together as "dry" deposition ( Table 1).

Inferential Method (IFM) Calculation procedure
The inferential method (IFM, also referred to as "concentration method", e.g., Peters and Eiden, 1992), calculates the dry deposition flux (DD; g m −2 s −1 ) as the product of the concentration (c; µg m −3 ) in the ambient air at a defined reference height (z) and a proportionality constant, the deposition velocity (vd; cm s −1 ) according to the following basic equation (Wesely and Hicks, 2000): Following Thimonier et al. (2019) we focused on NH 3 , NO 2 , HNO 3 (gases), and NH + 4 and NO − 3 (aerosols). As insitu measurements were only available for NH 3 and NO 2 , the calculation of DD based on vd and ambient air concentration was only applied for these substances, which usually account for a large part of TIN DD (Flechard et al., 2011). The contribution of HNO 3 , NH + 4 , and NO − 3 to DD was estimated following Schmitt et al. (2005). In detail, we proceeded as follows: (1) we obtained WD of NH + 4 and NO − 3 in the same way as for the CBM approach ( Table 1).
(2) We then calculated the DD of NH + 4, part , NO − 3, part , and HNO 3, gas based on an empirical relationship (Schmitt et al., 2005) from open field deposition, separately for broadleaved and coniferous forest. OD is not included here (Schmitt et al., 2005). (3) Gaseous DD of NH 3 and NO 2 was calculated by multiplying annual average ambient air concentrations with deposition velocities (see below). (4) As independent data to inform OD estimates was not available, we used the OD from the EBM approach (see below) also in the IFM approach ( Table 1). The contribution of this deposition pathway is usually very small. (5) The TIN TD estimate of the IFM is then calculated as the sum of the deposition fluxes from steps 1-4. (6) In order to allow a separate comparison of the deposition fluxes, we subtracted WD (from step 1) from the TD (step 5). Although this flux contains small parts of OD, we refer to it as dry deposition (DD) in subsequent steps of analyses (same for the other two methods).

Estimation of deposition velocities for NO 2 and NH 3
Previous reviews of deposition velocities for different forest types and substances reported a large variability of vd values, which appeared to be sensitive to several parameters such as receptor surface properties, meteorological conditions as well as seasonal, and diurnal variations (Hunová et al., 2016). To account for this variability, we derived two sets of deposition velocities. In the standard case, we used forest type specific vd values for each substance across Germany and the complete observation period. In a second case, we used semi-empirical site and season specific correction factors following Kirchner et al. (2014), to roughly account for the major spatial and temporal variations in vd.
In order to establish the fixed "forest type specific values" of the deposition velocity, we rely on four review studies, each covering several publications. We extracted the vd values for NH 3 and NO 2 reported to be best suited (median of vd's in one case) for the respective forest type by each of the four studies. We then aggregated these four values into one value per forest type and N species by using the midrange (average of the lowest and highest value). In order to roughly account for the range of uncertainty, we also calculate "standard value+30%" and "standard value-30%" ( Table 2).  The second set of deposition velocities was based on these "forest type specific values" but adjusted by five types of site specific correction factors, proposed for spruce stands in the Bavarian Alps (Kirchner et al., 2014) and an additional correction factor to account for the relevance of other tree species elsewhere in Germany.
vd cor = vd lit k sea k incl k wind k inv k up k tree (6) with: vd cor : corrected deposition velocity; vd lit : "forest type specific" deposition velocity based on literature references (cf. Table 2) and the correction factors: k sea : season; k incl : slope inclination; k wind : wind speed; k inv : inversion weather condition; k up : slope upwind and k tree : tree species. Dry deposition strongly depends on the season (Marner and Harrison, 2004;Mohan, 2016). Kirchner et al. (2014) therefore proposes the following seasonally dependent factors (k sea ): Spring (1.1); summer (1.2), autumn (1.0), and winter (0.8). A consideration of the different lengths of winter periods as proposed by Schmitt et al. (2005) has been omitted. The correction factor (k incl ) for the slope inclination (incl; %, see below for data source) is obtained according to Kirchner et al. (2014): k incl = 0.01 × incl+0.6. Parameterization of k wind is based on regionalized wind speed from 109 wind stations using a Leeward index according to Dietrich et al. (2019). Wind speed is known to have a strong effect on deposition fluxes (Lin et al., 1994;Gallagher et al., 1997;Erisman and Draaijers, 2003;Mohan, 2016). Based on the average wind speed in Germany at 10 m height of ∼3.4 m s −1 , the following classes with the corresponding correction factors were developed, assuming DD increases with wind speed due to a higher turbulence and transport into the forest: <1 m s −1 : 0.7; 1.0-1.9 m s −1 : 0.8; 2.0-2.9 m s −1 : 0.9; 3.0-3.9 m s −1 : 1.0; 4.0-4.9 m s −1 : 1.1; 5.0-5.9 m s −1 : 1.2; ≥6 m s −1 : 1.3. For meteorological inversions, Kirchner et al. (2014) proposed the following correction factors (k inv ) related to frequency of their occurrence: "rare": 1.0; frequent: 0.9; very frequent: 0.8. As no corresponding information is usually available for the investigated sites, we calculated the terrain exposure index (TEI) as predictor. The TEI is a terrain parameter, which indicates the degree to which a particular location is sheltered against advective flows and is thus particularly suitable to accumulation of cold air (Dietrich et al., 2019). The index values and correction factors are classified as follows: ≥1.2: 1.4; <1.2-≥1.1: 1.2; <1.1-≥1: 1.0; <1-≥0.9: 0.9 and <0.9: 0.8. Kirchner et al. (2014) have integrated upslope winds into their approach because they can periodically transport emissions from source regions (Benedict et al., 2013). However, the occurrence of slope upwinds and its influence on vd values (k up ) is difficult to parameterize. Due to the high surface roughness of forests, this effect was only considered for slopes >5 • . According to Kirchner et al. (2014), north-exposed slopes are less affected than south exposed ones. Accordingly, a k up factor of 1.1 was assigned for the celestial directions NW, N, NE, and E, for W and SE exposed slopes a factor of 1.2 and for SW and S exposed slopes a factor of 1.3. Vd values obtained from reviews are often stratified according to forest type (coniferous vs. deciduous). The approach by Kirchner et al. (2014) has been developed for spruce stands but only 34 of the 104 plots included in this study were pure spruce stands. The leaf area index (LAI) in spruce stands is often higher than in pine stands (Panferov et al., 2009;Goude et al., 2019). Zhang et al. (2003) found generally higher vd values for stands with higher LAI values. Corresponding, vd values for pine were assumed to be lower by a factor of 0.7 und for spruce to be 1.3 times higher compared to the "forest type specific values" for coniferous forest. For the same reasons (Bequet et al., 2011) we set a k tree of 0.9 for oak and 1.1 for beech trees. For all other trees k tree was set to 1.0.

Emission Based Deposition Model (EBM)
The approach for quantifying TIN deposition with the EBM is described in detail in Schaap et al. (2018). The deposition fluxes estimated by the EBM were provided by the German Environment Agency. In short, four major calculation steps are conducted in this model: (1) the chemical transport model LOTOS-EUROS (Schaap et al., 2008;Manders et al., 2017) is used to calculate DD as a product of modeled ambient air concentration fields of N species and modeled deposition velocities. Critical input data are meteorological data in high temporal resolution, spatial information of N emissions and receptor properties for dry deposition (e.g., land cover). (2) In the next step, modeled rain water concentrations from the LOTOS-EUROS model are used in combination with a few 100 stations of precipitation chemistry monitoring in Germany. These data serve to adjust the modeled rain water concentration distribution from the LOTOS-EUROS model using residual kriging. The generated concentration field is multiplied with high resolution precipitation data (1 × 1 km), to yield WD estimates ( Table 1).
(3) OD is calculated from fog water concentrations [estimated from previously determined rainwater concentration field and so called enrichment factors- Schaap et al. (2018)] in combination with cloud water deposition rates, which were calculated following the approach by Katata et al. (2008Katata et al. ( , 2011. In Katata et al. (2008) a simple linear regression for the fog deposition velocity has been derived from numerical experiments with a detailed multilayer land surface model (4) Finally WD, DD, and OD were combined in a spatial resolution of 1 × 1 km for the years 2000-2015. For each grid cell, land cover type specific deposition rates are available, including coniferous forest, deciduous forest, and mixed forest. Thus, deposition rates at the ICP Forests monitoring sites were extracted from the EBM results based on site coordinates and tree species / stand type (see Figure 1). In line with the CBM and the IFM approach, we refer to the sum of DD and OD estimates as "dry deposition."

Derivation of Large Scale Ambient Air Concentrations, Meteorological and Terrain Data
Differences in deposition estimates between the three approaches might be affected by a range of factors including meteorological and terrain characteristics as well as the level of ambient air concentrations. Daily data of temperature (T), solar radiation (RA), relative humidity (RH), wind speed (W), and precipitation (P) were obtained for the period from 1981 to 2015 using the observational data of the German Meteorological Service (Deutscher Wetterdienst, DWD). The regionalization of the daily meteorological data from the climate and precipitation stations of the DWD to the intensive monitoring plots was performed using the methods described in Dietrich et al. (2019). Terrain parameters [altitude, slope, aspect, terrain exposure index (TEI)] were also derived following Dietrich et al. (2019). The slope orientation of the plots was transformed into a continuous variable (aspect index) between 0 and 1 following Roberts and Cooper (1989).
Discrepancies between model estimates might also be affected by the degree of air pollution at the sites. To account for this aspect, we included the annual average NH 3 and NO X concentrations as predictors in the statistical models. Because in-situ measurements were only available at the much smaller subset of plots for which the IFM approach could be calculated, we relied on modeled ambient air concentrations to provide data for all plots over the entire observation period. In order to ensure independence from the chemical transport model used in the EBM approach (LOTOS-EUROS), we utilized data from the EMEP MSC-W model (Simpson et al., 2012). Annual average air concentrations based on emissions and meteorology for the years 2000-2015 according to the 2019 reporting status

Statistical Analysis of Spatial and Temporal Differences Between the Methods
We compared the three model types with summary statistics. For the IFM, we used the version with site specific correction factors vd cor . We also tested a version with site specific correction factors but excluding the tree species-specific correction (k tree ) and another version with the forest type specific vd only. To quantify the statistical association between the estimated deposition rates from the different methods, the root mean square error (RMSE), as well as the mean bias error (MBE), the mean absolute error (MAE), the coefficient of determination (R 2 ) and Legates and McCabe's efficiency (E 1 ) were used ( Table 3).
The E 1 (Legates and McCabe, 2013) provides additional information with E 1 = 1 indicating a perfect fit while E 1 = 0.0 indicates a model that is no better than the baseline comparison (the observed mean value-Null model).
Substantially flawed results are indicated by negative E 1 values (Legates and McCabe, 2013). We used a mixed effect model to relate the differences between deposition estimates ( TIN y,p ) from the three approaches to potential explanatory factors taking into account the "pseudoreplicated" structure of the data (same plots in different years) (Zuur et al., 2009). The model structure is as follows: were TIN y,p is the difference in TIN deposition in year y at plot p; b 0 : the intercept term; f 1 , f 2 ,...f n : unspecified (potentially) non-linear spline smoothing functions; x 1,yp , x 2,yp ,...,x n,yp : 1... n predictor variables in year y at plot p; Z p : a row in a model matrix including dummy variables for coding random effects for plots p, where p = 1,...,104; b p : a vector of random effects; ε: an independent and identically normally distributed error term with standard deviation e yp . Standard software to parameterize this type of model is available from the R (R Development Core Team, 2018) library mgvc (Wood, 2006), with additional calls to the libraries MASS (Venables and Ripley, 2003) and nlme (Pinheiro et al., 2018). The distribution of the response variable was checked with the function "fitdist" of the package "fitdistrplus v1.0-14" (Delignette-Muller and Dutang, 2015). In case of data deviating from normal distribution, logarithmic or square root transformations were performed with constants added where necessary to avoid negative values. All explanatory variables are summarized in Table 4. It should be noted that some of these variables (e.g., wind speed) are already included in the semi-empirical factors for IFM. Therefore, significant effects of these parameters could also indicate deficient parameterization of the IFM. We used high-resolution regionalized climate data instead of precipitation observations at the intensive monitoring plots as explanatory variables for the statistical models. This allows us to keep explanatory variables independent from input data of the three methods we aim to compare.
To identify the factors that influenced the differences in deposition estimates between the three approaches, we used a boosting framework called component-wise gradient boosting (Mayr et al., 2017), implemented in the R-package mboost (Hothorn et al., 2020). The preselection of potentially appropriate model variables takes into account parametric, non-parametric, spatial, and random effects (Bühlmann and Yu, 2003;Bühlmann and Hothorn, 2007). In mboost, the major tuning parameter of boosting is the number of iterations "mstop." To prevent overfitting, we used the implemented k-fold cross validation function to choose an appropriate number of boosting iterations. For the final model selection, we fitted the generalized additive mixed models (GAMM) with the preselected variables from the mboost procedure. Finally, all variables with a very low accumulated in-bag risk reduction (these values can be used to reflect variable importance) and nonsignificant (p < 0.05) effects were then stepwise removed from the model.

Magnitude and Variability of Total Inorganic Nitrogen Deposition
For the 1,237 annual observations considered in our study, the TIN TD averaged over time per plot ranges between 6.1 and 47.1 kg ha −1 a −1 and between 8.8 and 36.8 kg ha −1 a −1 with a mean of 20 and 18 kg ha −1 a −1 and a coefficient of variation of 36 and 26% for CBM and EBM, respectively ( Table 5). There was considerable variation in estimated DD for the CBM approach. The coefficient of variation calculated across all stands and years was 50%. As expected the data of the EBM approach shows a clearly lower variability (29%).

Comparison of the Canopy Budget and the Inferential Method
Since both methods use identical observations of WD, the comparison focuses on dry and total deposition. For both dry and total deposition, the use of forest type specific vd values without further corrections shows a large dispersion and discrepancy (E 1 , RMSE, MAE) between the two methods and the bias is relative high (Figure 2, Table 6). When site specific correction factors are taken into account, the agreement between the two methods considerably improves. The R 2 is relatively high, the E 1 rises substantially above zero, the RMSE and MAE are about halved, and also the MBE is low (−0.3 kg ha −1 a −1 ) for TD estimates ( Table 6).

Comparison of the Canopy Budget and the Emission Based Method
Estimates for the comparison of the CBM and EBM approaches are available for a large set of plots and years (Table 6, Figure 3). There is an overall good agreement (R 2 = 0.47) between CBM and EBM estimates for WD ( Figure 3A). On average, the wetonly corrected TIN bulk deposition (CBM) is 0.8 kg ha −1 a −1 lower than WD from EBM ( Table 6). For DD, large discrepancies between CBM and EBM appeared (R 2 = 0.02; E 1 = −0.05, Figure 3B). However, the bias of 2.8 kg N ha −1 a −1 is lower than for the EBM-IFM comparison (bias of 3.3 kg N ha −1 a −1 ). Due to the weak agreement of the DD estimates, the association between  Table 2). For vd min and vd max only the regression lines were displayed and not the individual points. All statistical parameters are given in Table 6.
the TD with R 2 = 0.16 is also very weak ( Figure 3C). The TIN TD calculated with EBM is lower on average by 2 kg N ha −1 a −1 compared with CBM. To compensate for interannual variability, also plot-specific mean values were compared ( Figure 3D). Using this aggregated values, the regression line approaches more closely the 1:1 line.

Comparison of the Inferential and the Emission Based Method
When comparing the wet-only corrected bulk open field deposition of the reduced dataset (i.e., only for those plots where IFM estimates are available) to WD of the EBM approach, the mean bias error (MBE) is slightly lower (−0.6 kg N ha −1 a −1 ) compared to the complete dataset ( Figure 4E and Table 6). In contrast, the DD from the IFM clearly shows higher values if calculated with site specific correction factors, compared to EBM (3.3 kg ha −1 a −1 , Figure 4B). As MBE of WD and DD compensate each other partly, the overall difference between the TIN TD estimates amounts to 2.7 kg N ha −1 a −1 . If IFM is calculated only with forest type specific deposition velocities, the association with the EBM approach is significantly deteriorated (Figure 4 and E 1 values in Table 6). For this version, the MBE increases to 6.9 and 6.3 kg N ha −1 a −1 for DD and TD, respectively.

Factors Influencing the Difference Between N Estimation Methods
The differences between the annual TIN fluxes calculated with the three methods were examined for possible influences of spatial, temporal, meteorological and relief factors. For each combination of models (CBM, IFM, EBM) and deposition pathway (WD, DD, TD) the best model as identified by the boosting procedure for variable selection is indicated in Table 7. The cross-comparison of the approaches is hindered by the unequal number of replicates. While 1,237 data points (deposition rates for specific plots-years) are available for the CBM-EBM comparison, only 194 are available for the IFM-EBM and IFM-CBM comparison, respectively. This is due to the limited number of available in-situ ambient air concentration measurements, which are required for the IFM approach. This means that differences in effects found for the CBM-EBM and the IFM-EBM comparisons, respectively, could be either caused by the difference in the amount and identity of input data or by the difference in methods (CBM, IFM). In order to evaluate this aspect, we repeated the analyses (GAMMs) for the CBM-EBM comparison, including only those observations, which are available for the IFM/EBM comparison (n = 194) ( Table 7). When comparing the results, the different spatial coverage must be taken into account. For example, there are no measurements of ambient NO 2 and NH 3 concentrations available from Bavaria. With the exception of slope inclination, orographic parameters do not actually appear in any of the models. There are several significant effects on TIN in the different models with varying effect strength. In the following, mainly effects with a high effect strength will be addressed, since they contribute most to the explained variance of the differences between the models. For example, in the CBM-EBM comparison the effect strength for precipitation is between 2 and −6 kg ha −1 a −1 (Figure 5A) whereas for the year of measurement this range is only about 1.5 kg ha −1 a −1 (Figure 5E, note the logarithmic scale). The TIN values for WD between CBM and EBM can only be explained to a small extent by the variables included in our analyses (R 2 = 0.36). The variables with the largest explanatory power are precipitation ( Figure 5A) and the location of the plots (Figure 6). The greater the amount of precipitation, either the EBM approach seems to systematically overestimate WD or the measurements of open field deposition used in the CBM and IFM are too low.
In addition, higher TIN WD rates by the EBM approach can also be observed in northwest Germany ( Figure 6A). In contrast, somewhat lower WD rates are estimated by the EBM compared to open field measurements in northeast and southwest Germany ( Figure 6A). Note that positive effects indicate a tendency for higher TIN deposition estimates of the CBM as compared to the EBM (Figures 5, 6).
For the comparison of CBM and EBM with respect to dry deposition estimates a larger number of influential variables were identified than for WD. With increasing temperature, wind speed and slope inclination, the difference of DD estimates between EBM and CBM tends to increase (Figure 5). A relatively high sensitivity on TIN is indicated for the partial effect of tree species (Figure 5F). While the differences are insignificant for all other tree species, they are significant for spruce forests ( Table 7). In contrast to WD, the spatial trends for DD are only weak (Figure 6B). Similarly, the temporal effect is only poor. However, especially since 2012, there has been a slight downward trend, i.e., a tendency for higher values of the EBM compared to the CBM ( Figure 5E).
For the comparison of the TIN total deposition between the IFM and the CBM (possible for 194 plot-years), the variable selection algorithm suggested the location of plots and tree species as the sole potentially relevant explanatory variables. However, these effects turned out to be insignificant in the GAMM, with a corresponding low explanatory power (R 2 = 0.12).
As for the CBM-IFM comparison, the IFM-EBM comparison was also limited to the 194 plot-years with observed ambient air quality (Figure 7). For WD, the statistical model identified year of measurement and a spatial effect as explanatory variables. In contrast to the comparison of CBM and EBM a significant effect of precipitation rates was not found. It should be noted that there are only very few stations with very high annual precipitation rates in the IFM subset of plots, e.g., the Bavarian Alps are missing. Although a precipitation effect on differences  Table 6. in WD was invisible, the statistical model identified a tendency for higher TD estimates of the EBM compared to IFM for plotyears with high regionalized precipitation rates. Another effect identified for the IFM-EBM comparison is a systematically lower TIN DD estimate of EBM compared to the IFM for spruce sites (Figure 7E). This effect of tree species is in agreement with the results of the larger set of plots used for the EBM-CBM comparison (see above). The IFM-EBM comparison also indicates an effect of solar radiation, NO x and NH 3 air concentrations (based on data from the EMEP model) for dry deposition. However, the prognosis intervals show a high uncertainty. The comparison of the IFM and EBM methods does not reveal any spatial pattern for TD, which is likely due to the limited sample size.  Table 2). For vd min and vd max only the regression lines were displayed and not the individual points. All statistical parameters are given in Table 6. the respective variants ( Table 7). For example solar radiation is significant for IFM-EBM but not for the other variants. However, this can also be an effect of the limited sample size. It can be seen for TD, that when comparing CBM and EBM (n = 1,237), all variables listed in the table are significant, with the exception of solar radiation. -Spatial and temporal effects are mainly shown for the comparison between CBM and EBM. This is likely also the result of the much larger sample size. Spatial effects with a large effect strength are mainly apparent for WD. -With the exception of slope inclination, orographic parameters do not actually appear in any of the models.

DISCUSSION
We compared three commonly used models for the estimation of TIN deposition to forests: a canopy budget model (CBM, Ulrich, 1994), the inferential method (IFM) based on observations of ambient air concentrations and an emission-based-approach (EBM) using a chemical transport model (LOTOS-EUROS).
The relatively large number of plot-years included in this study allowed us to identify general patterns of differences between the three approaches. At the same time, the amount of data made it impossible to analyze each plot-year in detail. Therefore, although several important explanatory variables were included in the statistical models, some patterns occurring at specific subsets of the data might not be well-represented in the results. The discussion comprises two main sections: (1) uncertainties in methods and measurements; (2) comparison of TIN deposition estimates between the three methods.

Uncertainties in Methods and Measurements
The three methods to estimate TIN deposition compared in this study (CBM, IFM, and EBM) are frequently applied because they can be used with relatively low effort compared to methodologically advanced micrometeorological methods. This also implies, however, that each method suffers from substantial limitations. Ellermann et al. (2018) estimated the uncertainty on nitrogen deposition for Danish land areas to be 27-43%. They interpret the high uncertainty as a result of partial uncertainties of the various N species that contribute to TIN TD. In this study, the uncertainty of the different methods can only be quantified very roughly, as the results are obtained from a combination of measurements and models with different input data and many assumptions. The following aspects of uncertainty for the three models should be considered when interpreting the results. For WD measurements alone, uncertainties of 10-40% can be assumed due to errors in sampling and analysis (Zimmermann et al., 2006). The higher small-scale variability in forest stands is taken into account by a higher number of collectors, therefore the uncertainty in throughfall deposition is expected to be of comparable, but potentially somewhat larger, magnitude. A further uncertainty results from the application of uniform conversion factors to convert from bulk to wet deposition. These factors do not only affect wet deposition, but also the TABLE 7 | Significance and direction of the effects explaining differences between TIN deposition estimates for different deposition pathways (DP) and pairs of models. R 2 , adjusted coefficient of determination (R 2 ); N, number of site-years. Asterisks indicate significance levels (*p < 0.05, **p < 0.01, ***p < 0.001, n.s., not significant but included in model, see Materials and Methods). Blank fields indicate that the parameter has not been considered relevant by the variable selection algorithm (mboost) or that it was not significant in final mixed effect models. Parameters that were not included in any model are not listed. TIN model = Combination of models for which differences in TIN deposition estimates are analyzed; C = coordinates (easting and northing, Gauss-Krüger); Y: year of deposition data (measurement or estimated); TR = tree species (sp = sig. differences for spruce; pi = sig. differences for pine; mi= sig. difference for mixed stands). Effects must be interpreted relative to deciduous forest; P, annual precipitation sum, [mm]; T, annual mean temperature, [ • C]; RH, annual mean relative humidity, [%]; RA, annual solar radiation sum, [MJ m −2 ]; W, annual mean wind speed, [m s −1 ]; S, slope inclination, [ • ]; NO X , NO X concentration, [µg N m −3 ]; NH 3 , NH 3 concentration, [µg N m −3 ]; VT, variant without tree species specific correction; ↑, increasing (not necessarily linear) difference between models with increasing parameter value; ↓, decreasing (not necessarily linear) difference between models with increasing parameter value; ∼, multi-directional effect (e.g., first increasing, then decreasing).

DP N TIN models C Y TR P T RH RA W S NO
results of the CBM method and the estimation of particulate deposition from observed wet deposition for the IFM method. The national average values we used (0.95 for NH + 4 , 0.9 for NO − 3 ) have a standard deviation of 0.25 and 0.22 for NH + 4 and NO − 3 , respectively. They are based on 79 (NH + 4 ) and 86 (NO − 3 ) plot-years of parallel wet and bulk sampling across Germany (Gauger et al., 2008). Bulk to wet conversion factors can vary greatly between regions and are particularly uncertain for ammonium, where losses from bulk collectors have been reported (Stedman et al., 1990;Fürst, 2016). The canopy budget model is also affected by the wet-only conversion for sodium (tracer substance in the CBM). We used the national average of 0.81 for Na + , but the spatio-temporal variability for this value across Germany is very high (standard deviation = 0.2; Gauger et al., 2008). Further studies are required to identify variables that could explain the spatio-temporal variability and could accordingly allow for a regionalization of the conversion factors.
The CBM approach according to Ulrich (1994) relies on the robust measurements of six flux rates (Na + , NH + 4 , NO − 3 for both open field and stand precipitation). While other canopy budget models (Draaijers and Erisman, 1995;De Vries et al., 2001) require even more parameters and assumptions, the accurate monitoring of these fluxes is still challenging. Contamination with biological material can alter the composition especially of throughfall and stemflow, which is sometime challenging to detect despite the rigorous QA/QC rules applied (Mosello et al., 2005). Furthermore, the CBM approach relies on the assumption of similar concentration ratios (NH + 4 :Na + and NO − 3 :Na + ) in interception deposition and wet deposition. It has primarily been formulated for areas where cloud water droplets dominantly contribute to interception deposition and where substantial amounts of the airborne tracer substance (Na + ) are present (Ulrich, 1994). Deviations from these conditions are expected for the majority of the German intensive monitoring plots. The resulting uncertainties in deposition estimates remain to be investigated. Comparisons between micrometeorological methods and CBM have reported high deviations (Gallagher et al., 1997;Mohr et al., 2005). Another weakness of the CBM approach according to Ulrich (1994) is that it conceptually underestimates TIN DD as net uptake of TIN in the canopy compensates DD fluxes. Uptake of TIN in the canopy may occur via stomatal uptake of gaseous N species by the trees or via ion exchange. The conversion of inorganic to organic N species by microorganisms in the canopy additionally contributes to the underestimation of TIN DD. Canopy budget models with explicit parameterization of canopy uptake (Draaijers and Erisman, 1995;De Vries et al., 2001) estimate on average higher TIN deposition rates (Mohr et al., 2005). Flechard et al. (2011) demonstrated that discrepancies between N deposition estimates from different IFM models at one site can be several times larger than between sites. One possible reason for this is that the IFM approach usually is a combination of observations and modeling with numerous assumptions. This also applies to the IFM approach used in our study. The gaseous components NO 2 and NH 3 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 parametrization of deposition velocity vd. Several mechanistic models for vd exist, which e.g., utilize detailed information on meteorology and soil water availability to model stomatal opening. This allows for the parametrization 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 dry deposition fluxes, additional uncertainty arises from assumptions about input parameters. We chose a different approach, starting from typical deposition velocities from the literature ("forest type specific values"). These forest type specific values were adjusted using the most important determinants for vd in order to derive "semi-empirical correction factors." The results show a much better agreement with the two other methods (CBM, EBM) than the uncorrected version. In order to account for other tree species at a majority of the intensive forest monitoring sites in Germany, we extended the IFM approach established by Kirchner et al. (2014) at 9 spruce sites in the Bavarian Alps by the tree species specific correction term. This extension and the use of the terrain exposure index (TEI) as a proxy for the frequency of inversion weather conditions were the only modifications made. Applying these semi-empirical correction factors reduced the bias between the IFM and both the EBM and the CBM approach substantially ( Table 6). Additionally, in order to illustrate the uncertainty caused by the forest type specific deposition velocity values, the results for 30% higher and lower vd, which roughly corresponds to the variability of vd from the literature, were integrated as regression lines (Figures 2, 4).
Since the deposition of gaseous HNO 3 and particulate NO − 3 and NH + 4 often only accounts for a small proportion of dry deposition (Thimonier et al., 2019) and no observations are available, we used an empirical relationship (Schmitt et al., 2005) which can be parameterized by forest type and bulk deposition of TIN. The equations were obtained from 77 monthly measured deposition values of the different N species and the coefficient of determination was 0.92. We performed a limited test whether the approach can be transferred to other deposition and meteorological regimes based on the deposition data reported by Thimonier et al. (2019), which resulted in an underestimation of about 1 kg N ha −1 a −1 . However, the test is only valid to a limited extent. Data used for model parameterization and evaluation come from different years but partly from the same sites. Furthermore, only bulk deposition was measured at these sites, so the average conversion factors for Germany to estimate WD (Gauger et al., 2008) had to be used. Therefore, we can only assume that a rough estimate of gaseous HNO 3 and particular NO − 3 and NH + 4 deposition is possible with this approach.
The EBM approach contains uncertainties at very different methodological and spatial levels. Comparison of spatial data in a 1 * 1 km grid resolution as for the EBM with point data as for the CBM and IFM should be done very carefully, especially in areas with complex orography or in regions with a fine-grained land cover pattern (edge effects) and with a high variability of N emissions. Actual deposition rates at a specific location may differ from the average of the corresponding EBM grid cell, although EBM output is reported on a land cover specific level (conifer, broadleaved and mixed forest). Similar to the other two approaches, EBM is affected by the uncertainties in wet deposition assessments. More importantly, as precipitation is highly variable in space and time, further uncertainty originates from the regionalization of precipitation, depending on regionalization strategy, orographic conditions and scale. Considering these challenges to correctly estimate TIN WD, the average deviation of −0.8 kg N ha −1 a −1 compared to wet-onlycorrected bulk deposition measurements for 1,237 plot-years in Germany is remarkably low. For dry deposition estimation Saylor et al. (2019) pointed out that the algorithms used in atmospheric chemistry models to predict particle deposition velocity are highly uncertain. In particular, estimates for forests show a weak agreement with available measurements (Saylor et al., 2019). Furthermore, the performance of EBM depends on the spatial resolution and quality of their input data, in particular the emission inventory. Meteorological data are available in a high temporal, but low spatial resolution of 7 × 14 km . Accordingly, Simpson et al. (2011) estimated an error of TIN deposition estimates of 30% for the different regional models and approaches in Europe. A conceptual problem faced by all types of models is that relations established at specific locations are extrapolated to other areas (e.g., assumption of the CBM, dry deposition velocities, scavenging ratios, fog water enrichment factors) and proper parametrization of the transfer of these relations to other locations is challenging.

Comparison of TIN Deposition Estimates Between the Three Methods
We found on average a difference of −0.3 kg ha −1 a −1 and +3.3 kg ha −1 a −1 between the CBM and IFM approaches, with and without applying the extended semi-empirical deposition velocity correction factors given by Kirchner et al. (2014), respectively. The correction improves the fit between CBM and IFM across Germany (lowlands and mountain ranges), although it has been developed for the Bavarian Alps. Schmitt et al. (2005) found slightly higher estimates of the IFM (median difference +2.4 kg N ha −1 a −1 ) when compared to a different CBM (De Vries et al., 2001), which, however, usually yields higher estimates of TIN deposition, compared to the Ulrich (1994) model used in our study (Mohr et al., 2005). Zimmermann et al. (2006) found a good agreement between the two methods, while Kirchner et al. (2014) stated that the two methods deviated at all sites. A slightly lower TIN deposition estimated by the CBM compared to the IFM may be expected due to the conceptual underestimation of TIN deposition by the specific type of CBM we used (Ulrich, 1994). However, the regression functions from Schmitt et al. (2005) used in the IFM approach for estimating HNO 3 , NH + 4 , and NO − 3 deposition may also lead to additional uncertainty of the deposition rates.
The statistical analysis of the differences in estimated TIN deposition between the IFM and CBM approaches found no significant effects of tree species, measurement year, coordinates, altitude, slope, aspect, aspect index, terrain exposure index (TEI), wind speed, air concentrations or precipitation. A large effect of the random factor "site" in the mixed model might indicate that further site-specific properties exists which contribute to the observed differences between model estimates in addition to random errors. Stand height, canopy closure, frontal area index (total area of roughness elements projected in the wind direction per unit ground area), crown density, and/or stand volume could be some of these site specific factors affecting local deposition processes (Erisman and Draaijers, 2003;Vesala et al., 2005;De Schrijver et al., 2007;Nakai et al., 2008).
A few studies exist where CBM and/or IFM have been compared to micro-meteorological measurements of higher accuracy. Recently, Brümmer et al. (2020) have tested an eddy covariance approach for total reactive nitrogen in a forested area of the Bavarian forest with low levels of N deposition. For the period 2016-2018 they report a measured TIN DD rate of around 4.4 kg ha −1 a −1 compared to estimates of 5.2 and 6.9 kg ha −1 a −1 from the same EBM as used in our study (based on an uncorrected and an improved classification of land use types, respectively). In parallel, a CBM approach is conducted at the site since several years (Beudert and Breit, 2014). However, the CBM results are not comparable with our study, because the variants of the CBMs they use clearly differ from the Ulrich (1994) model used in our study. Mohr et al. (2005) found clearly lower N fluxes (−27 kg N ha −1 a −1 ) with the CBM compared to a micrometeorological method at a site exposed to high ambient air levels of NH 3 .
The comparison of WD using the EBM approach with the CBM shows a fairly good agreement. However, there seems to be a systematic deviation for estimates with increasing precipitation rates ( Figure 5A). This effect mainly originates from two mountainous stations. At these stations, the simulated precipitation amount clearly exceeds the open field precipitation observations. If the corresponding stations are excluded from the statistical analysis, the effect is still significant, but the effect strength is somewhat smaller (results not shown). The deviation in the precipitation data might be caused by a systematic underestimation of precipitation amount for observations with a large proportion of snow or strong rain events (sampler overflow) and/or low performance of regionalization models for precipitation in orographic complex landscapes as used for WD estimation (1 * 1 km). The elevated deviance between EBM and CBM estimates in some regions might be attributable to the uniform application of bulk to wet conversion factors ( Figure 6A).
The EBM approach estimates on average lower DD compared to CBM and IFM approaches ( Table 6). This pattern is more distinct for spruce plots (Figure 5F, Figure 7E) and does not diminish if the tree-species correction factor of the IFM approach is deactivated (Table 7, model "IFM VT -EBM"). As the CBM and the IFM rely on independent approaches to estimate the dry deposition, this might indicate that the EBM deposition estimates are less reliable for spruce plots. When interpreting the partial effects, however, the relationship between individual variables must always be taken into account. Accordingly, the change in nitrogen deposition should not be deduced directly from the tree species spruce, for example, but rather the respective precipitation conditions (and other variables) at the site must also be considered. Several meteorological parameters significantly affect the differences between DD estimates, which are also known to have an effect on deposition rates [e.g., via stomatal resistance, particle deposition velocity; Han et al. (2011), Mohan (2016)]. On the other hand, the different terrain parameters hardly contribute to an explanation. This may be due to incorporation of terrain information into the regionalization of the climatic variables that were used as predictors (Dietrich et al., 2019).
In a previous comparison between the EBM used in our study and the CBM according to De Vries et al. (2001), TIN TD estimates were reported to roughly agree . At some locations, however, differences of up to 40% were found . Our analysis is, however, not directly comparable to the Schaap et al. (2018) study, as (1) a different type of CBM was used; (2) a larger number of plots with an extended regional coverage of Germany and longer observation periods are included; (3) a bulk-to-wet correction of NH + 4 , NO − 3 , and Na + fluxes has been applied prior to CBM calculations as recommended by Adriaenssens et al. (2013) in order to harmonize the WD calculation between the three approaches. This correction increases the estimated TIN TD of the CBM by ∼1.2 kg ha −1 a −1 on average.

CONCLUSIONS AND OUTLOOK
The demand for N deposition estimates with high accuracy as well as large spatial and temporal resolution is growing. Forests pose a special challenge in this regard, due to the high importance of dry deposition, which is more difficult to quantify. In the framework of long-term monitoring at ICP Forests sites in Germany, we compared three methods for the estimation of total inorganic nitrogen deposition to forests (CBM, IFM, EBM). If all approaches would yield accurate results, we would have expected similar TIN deposition estimates from the IFM and the EBM and lower values from the CBM, as it includes a conceptual underestimation. Contrarily, we found the EBM provided on average lowest estimates. The deviation between the EBM and both of the other methods was especially pronounced for the dry deposition at spruce plots. Differences of dry deposition estimates between all methods were found to be affected by meteorological conditions, which are known to regulate deposition velocity. The average discrepancy in TIN deposition according to the method yielding on average highest deposition rates (IFM) and the method suggesting on average lowest deposition rates (EBM) was 6.3 kg N ha −1 a −1 (uncorrected IFM) or 2.7 kg N ha −1 a −1 (IFM with site specific corrections). We hypothesize that a combination of different factors contribute to the discrepancies between the three approaches. They include the comparison of observations at the plot scale (CBM, IFM) with gridded information (EBM, based on a 1 × 1 km resolution and coarser for some calculation steps), conceptual limitations of the utilized CBM version, and varying accuracy of the EBM for different forest types.
As all methods are subject to considerable uncertainty, we cannot conclude whether any of the methods provides more or less accurate estimates. Instead, the results can only indicate aspects worthwhile to consider for future methodological improvements. Recent developments contributing to a more accurate quantification of TIN deposition to forests include for example (i) continuous improvements of EBM and the underlying emission inventories (Schaap et al., 2015(Schaap et al., , 2017, (ii) networks of low-cost monitoring stations for ambient air concentrations and meteorological conditions improving data availability (e.g., Karagulian et al., 2019;Weissert et al., 2020), and (iii) ongoing research in CBM and similar approaches, like surface-wash and surrogate surface sampling (Aguillaume et al., 2017;Karlsson et al., 2019). In addition, more high-accuracy validation datasets (e.g., Brümmer et al., 2020) are required to quantify the performance of the different modeling approaches.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request from the institutions holding the data: Climate and terrain data by Institute of Geography, University of Hamburg, Germany. Emission based deposition data by UBA, Dessau, Germany. Measured deposition and air concentrations by Thünen Institute of Forest Ecosystems, Eberswalde, Germany. Deposition calculated by inferential methods by Northwest German Forest Research Institute (NW-FVA). Requests to access these datasets should be directed to Jan Wehberg, jan.wehberg@uni-hamburg.de; Markus Geupel, markus.geupel@uba.de; Anne-Kartin Prescher, anne.prescher@thuenen.de; Bernd Ahrends, bernd.ahrends@ nw-fva.de.

AUTHOR CONTRIBUTIONS
BA and AS: designed the study, developing of methods, preparation of deposition data, and data analysis. JW: regionalization of climate data and realization of terrain analysis. A-KP: preprocessing and aggregation of air concentration data. MG: emission based deposition data. HA and HM: supervising the work. All authors contributed to writing.

FUNDING
The evaluation was based on data that was collected by partners of the official UNECE ICP Forests Network (http://icpforests.net/contributors). 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 Systeme (FutMon).