Satellite-derived leaf area index and roughness length information for surface–atmosphere exchange modelling: a case study for reactive nitrogen deposition in north-western Europe using LOTOS-EUROS v2.0

The nitrogen cycle has been continuously disrupted by human activity over the past century, resulting in almost a tripling of the total reactive nitrogen fixation in Europe. Consequently, excessive amounts of reactive nitrogen (Nr) have manifested in the environment, leading to a cascade of adverse effects, such as acidification and eutrophication of terrestrial and aquatic ecosystems, and particulate matter formation. Chemistry transport models (CTMs) are frequently used as tools to simulate the complex chain of processes that determine atmospheric Nr flows. In these models, the parameterization of the atmosphere–biosphere exchange of Nr is largely based on few surface exchange measurement and is therefore known to be highly uncertain. In addition to this, the input parameters that are used here are often fixed values, only linked to specific land use classes. In an attempt to improve this, a combination of multiple satellite products is used to derive updated, time-variant leaf area index (LAI) and roughness length (z0) input maps. As LAI, we use the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD15A2H product. The monthly z0 input maps presented in this paper are a function of satellitederived normalized difference vegetation index (NDVI) values (MYD13A3 product) for short vegetation types (such as grass and arable land) and a combination of satellite-derived forest canopy height and LAI for forests. The use of these growth-dependent satellite products allows us to represent the growing season more realistically. For urban areas, the z0 values are updated, too, and linked to a population density map. The approach to derive these dynamic z0 estimates can be linked to any land use map and is as such transferable to other models. We evaluated the sensitivity of the modelled Nr deposition fields in LOng Term Ozone Simulation – EURopean Operational Smog (LOTOS-EUROS) v2.0 to the abovementioned changes in LAI and z0 inputs, focusing on Germany, the Netherlands and Belgium. We computed z0 values from FLUXNET sites and compared these to the default and updated z0 values in LOTOS-EUROS. The root mean square difference (RMSD) for both short vegetation and forest sites improved. Comparing all sites, the RMSD decreased from 0.76 (default z0) to 0.60 (updated z0). The implementation of these updated LAI and z0 input maps led to local changes in the total Nr deposition of up to ∼ 30 % and a general shift from wet to dry deposition. The most distinct changes are observed in land-use-specific deposition fluxes. These fluxes may show relatively large deviations, locally affecting estimated critical load exceedances for specific natural ecosystems. Published by Copernicus Publications on behalf of the European Geosciences Union. 2452 S. C. van der Graaf et al.: Satellite-derived leaf area index and roughness length for nitrogen deposition


Introduction
The nitrogen (N) cycle has been continuously disrupted by human activity over the past century (Fowler et al., 2015;Galloway et al., 2004;Galloway et al., 2008), resulting in a doubling of the total reactive nitrogen (N r ) fixation globally and even a tripling in Europe. As a result, excessive amounts of N r , defined as all N compounds except N 2 , have manifested in the environment contributing to acidification and eutrophication of sensitive terrestrial and 40 aquatic ecosystems (Bobbink et al., 2010a;Paerl et al., 2014). NO x and NH 3 affect air quality through their significant role in the formation of particulate matter, impacting human health and life expectancy (Lelieveld et al., 2015;Bauer et al., 2016;Erisman and Schaap, 2004). N r also influences climate change through its impact on greenhouse gas emissions and radiative forcing Butterbach-Bahl et al., 2011). As N r forms are linked through chemical and biological conversion in one another within the environmental compartments, one atom 45 of N may even take part in a cascade of N r forms and effects (Galloway et al., 2003). To minimize these adverse effects, effective nitrogen management and policy development, therefore, require consideration of all N r forms simultaneously.
With the scarceness and inadequate distribution of available ground measurements, especially for reduced N r , the most important method to assess and quantify total N r budgets on a larger spatial scale to date remains the use of 50 models. Models, chemistry transport models, in particular, are used for understanding the atmospheric transport and the atmosphere-biosphere exchange of nitrogen compounds. Most chemistry transport models compare reasonably with observations for oxidized forms of N r , but need improvement when it comes to the reduced forms of N r (Colette et al., 2017). Modelled NH 3 fields are in general uncertain due to the highly reactive nature and the uncertain lifetime of NH 3 in the atmosphere. More importantly, NH 3 emissions that are used as model input are very 55 complex to estimate and remain highly uncertain (Reis et al., 2009;Behera et al., 2013), for example, due to the diversity in NH 3 volatilization rates originating from different agricultural practices. Recently, a lot of effort has been made to improve the spatiotemporal distributions of bottom-up NH 3 emissions (e.g. Skjøth et al., 2011)). Emissions can also be estimated top-down through the usage of data assimilation and inversion techniques. Optimally combining observations and chemistry transport models have already enabled us to 60 create large-scale emission estimates for various pollutants (e.g. (Curier et al., 2014;Abida et al., 2017)), for instance for NO 2 , and will likely also be used for large-scale NH 3 emission estimates in the future.
Most data assimilation and inversion methods rely on the assumption that sink terms in the model hold a negligible uncertainty. To obtain reasonable top-down emission estimates, we must thus also aim to reduce the uncertainty involved on this side of models. The sink strengths of trace gases and particles in chemistry transport models are realistic, spatial-and time-variant LAI and z 0 values that are derived from optical remote sensors.
The LAI is defined as the one-sided green leaf area per unit surface area (Watson, 1947). The LAI serves as a measure for the amount of plant canopy, and herewith directly related to energy and mass exchange processes. As a result, the LAI is nowadays used as one of the main parameters in many ecological models. In deposition modelling, stomatal uptake is often parameterised using the LAI. The LAI can be determined in the field using either direct 80 methods, such as leaf traps, or indirect methods, such as hemispherical photography (Jonckheere et al., 2004). For larger areas, the LAI can be simulated using land surface or biosphere models. Another group of indirect methods to estimate the LAI for large regions is the use of optical remote sensing. The LAI can, for instance, be estimated using empirical relationships between LAI and vegetation indices (e.g. (Soudani et al., 2006;Davi et al., 2006;Turner et al., 1999)) or by inversion of canopy reflectance models (e.g. (Houborg and Boegh, 2008;Myneni et al., 2015)). A well-85 known example of the latter is the LAI product from the MODIS instrument, which we will use in this study.
The z 0 is used to describe the surface roughness. The surface roughness serves as a momentum sink for atmospheric flow and is, therefore, an important term in atmospheric modelling. The interaction between the boundary layer and the roughness of the Earth's surface results in shear stress that affects the wind speed profile. Under neutral conditions the resulting wind profile can be approximated using a logarithmic profile: 90 ( ) = * (eq. 1) ( ) represents the mean wind speed, * the friction velocity and the Von Kármán constant. Here, is a constant that represents the height at which the wind speed theoretically becomes zero. The z 0 can be estimated from in-situ wind speed measurements using bulk transfer equations. More recently, several studies have shown that z 0 for specific, uniform land cover types can also be estimated from optical remote sensing measurements, for instance 95 using vegetation indices (e.g. (Xing et al., 2017;Yu et al., 2016;Bolle and Streckenbach, 1993;Hatfield, 1988;Moran, 1990)). The z 0 can also be estimated using (satellite-derived) vegetation height (e.g. (Raupach, 1994;Plate, 1971;Brutsaert, 2013;Schaudt and Dickinson, 2000)).
The use of optical remote sensing data holds promising potential for improvements of the representativeness of the surface characterization in chemistry transport models. Here, we present an approach to derive monthly z 0 input maps using satellite-derived NDVI values (MYD13A3 product) for short vegetation types and a combination of satellite-derived forest canopy height and LAI for forests. We validate these z 0 values by comparing them to z 0 values computed from FLUXNET observations. We also update the z 0 values for urban areas, using a population density map. We use the updated z 0 values, as well as the MODIS-LAI, as input in LOTOS-EUROS to illustrate the effect on transport and deposition modelling of N r components. We evaluate the sensitivity of the N r deposition 105 fields to these input parameters, focusing on Germany, the Netherlands and Belgium. Moreover, we quantify and present the implications for land use specific fluxes on a model subpixel level. Also, we compare our model outputs with wet deposition measurements of NH 4 + and NO 3 and surface concentration measurements of NH 3 and NO 2 .

Model description
The LOTOS-EUROS model is a Eulerian chemistry transport model that simulates air pollution in the lower troposphere (Manders et al., 2017). In this study the horizontal resolution is set to 0.125⁰ by 0.0625⁰, corresponding to pixels of approximately 7 by 7 kilometres in size. The model uses a five-layer vertical grid that extends up 5 km above sea level, starting with a surface layer with a fixed height of 25 meters. The next layer is a mixing layer, 115 followed by two time-varying dynamic reservoir layers of equal thickness, and a top layer up to 5 km. LOTOS-EUROS follows the mixed layer approach and performs hourly results using ECMWF meteorology (European Centre for Medium-Range Weather Forecasts, 2016). The gas-phase chemistry uses the TNO CBM-IV scheme (Schaap et al., 2009) and the anthropogenic emissions from the TNO-MACC-III emission database (Kuenen et al., 2014). The wet deposition parameterisation is based on the CAMx approach, and includes both in-cloud and below-120 cloud scavenging (Banzhaf et al., 2012). LOTOS-EUROS makes use of the CORINE/Smiatek land use map to determine input values for surface variables.

Dry deposition
The dry deposition flux of gases is computed following the resistance approach, in which the exchange velocity is equal to the reciprocal sum of the aerodynamic resistance , the quasi-laminar boundary layer resistance and 125 the canopy resistance : = (eq. 2) and are both influenced by the wind profile, which is computed with (eq. 1). The wind profile, in turn, depends on roughness length z 0 associated with different land use classes. The aerodynamic resistance is computed as follows: Here, ℎ is the canopy height, which is pre-defined per land use class. The empirical function Φ is taken from Businger et al. (1971) and depends on the state of the atmosphere. Depending on the value of Monin-Obukhov length , the following equations are used: For a stable atmosphere ( > 0): For an unstable atmosphere ( < 0): For a neutral atmosphere is it equal to unity. The Monin-Obukhov length is dependent on z 0 , and is determined as follows: Here, , , , and are constants (0.004349, 0.003724, -0.5034, 0.2310 and -0.0325 respectively) and = −0.5 (3.0 − 0.5 | |) with near-surface wind speed and exposure factor , depending on cloud cover and solar zenith angle. The wind speed at a reference height (10 meters) is used to obtain the friction velocity * : * = / (eq. 7) The quasi-laminar boundary layer resistance is a function of the cross-wind leaf dimension and the wind speed at canopy top, (ℎ), following the parameterisation presented in McNaughton and Van Den Hurk, 1995: is set to 0.02 m for arable land and permanent crops, and to 0.04 m for deciduous and coniferous forests. For other land use classes , and subsequently, the , is equal to zero. The canopy resistance is computed using the DEPAC3.11 (Deposition of Acidifying Compounds) module (van Zanten et al., 2010). is a parallel system of the resistances of three different pathways, the external leaf surface or cuticular resistance , the effective soil 150 resistance , and the stomatal resistance , and is defined as: The external leaf surface resistance is a function of the surface area index (SAI) and the relative humidity. The SAI is a function of the LAI. The effective soil resistance , is the sum of the in-canopy resistance and the soil resistance . Soil resistance has a fixed value, depending on the land use class and conditions 155 (frozen, wet or dry). In case of * > 0 the in-canopy resistance for arable land, permanent crops and forest is computed as follows: = * (eq. 10) For * < 0 and other land use classes a fixed value is used. The stomatal conductance for optimal conditions is the product of the LAI and the maximum leaf resistance from (Emberson et al., 2000). This maximum value is reduced 160 by correction factors for photoactive radiation, temperature and vapor pressure deficit to obtain the stomatal conductance . The is then equal to 1/ . The resistance parameterizations differ with land us type. A total of nine different land use classes are used in DEPAC. LOTOS-EUROS uses a fixed z 0 value for each of these land use classes. The default LAI values are also linked to the DEPAC classes, and follow a fixed temporal behaviour that describes the growing season of each land use class (Emberson et al., 2000). The bi-directional exchange of NH 3 is

Datasets
The following section gives a short description of all the datasets that are used in this study. Firstly, a description of the LAI dataset is given. Subsequently, the datasets that are used to derive the updated z 0 maps are described.
Finally, the in-situ observations used for validation are discussed in the last paragraphs.

MCD15A2H Leaf Area Index
The satellite-derived Leaf Area Index (LAI) is a combined product of the MODIS instruments on board the Terra and Aqua satellites (Myneni et al., 2015). The LAI algorithm compares bidirectional spectral reflectances observed 175 by MODIS to values evaluated with a vegetation canopy radiative transfer model that are stored in a look-up table.
The algorithm then archives the mean and the standard deviation of the derived LAI distribution functions. We used the 6 th version of the product, MCD15A2H, which has a temporal resolution of 8 days and a spatial resolution of 500 meters.

180
The Normalized Difference Vegetation Index (NDVI) is a vegetation index computed with reflectances observed by the MODIS instrument on board of the Aqua satellite (Didan, 2015). The NDVI is the normalized transform of the near infrared to the red reflectance and is expressed by: We used the MYD13A3 product, which is the monthly NDVI product with a spatial resolution of 1 km.

Forest canopy height
The forest canopy height is derived from LIDAR (Light Detection And Ranging) data acquired by the GLAS (Geoscience Laser Altimeter System) instrument aboard the ICESat (Ice, Cloud, and land Elevation Satellite) satellite (Zwally et al., 2002). This instrument was an altimeter that transmitted a light pulse of 1024 nm and recorded the reflected waveform. We used the global forest canopy height product developed by Simard et al.
(2011), which has a spatial resolution of 1 km.

Population density map
The population density grid used in this study available for all European countries and provided by the European Environmental Agency (Gallego, 2010). The population density is disaggregated with the CORINE Land Cover inventory of 2000, using data on population per commune. The resulting population density grid is downscaled to a 195 spatial resolution of 100 meters.

CORINE Land Cover
The CORINE Land Cover (CLC) is a land use inventory that consists of 44 classes (European Environmental Agency, 2014). The CLC is produced by computer-assisted visual interpretation of a collection of high resolution satellite images. The CLC has a minimum mapping unit of 25 ha and a thematic accuracy of >85%. We used the 200 latest version of the product, CLC 2012, in this study.

In-situ reactive nitrogen observations
The modelled NH 4 + and NO 3 wet deposition fluxes are compared to observations of wet-only samplers. We used observation from the Dutch LMRe network (Van Zanten et al., 2017) and the German Lander network (Schaap et al., 2017).The location of the stations can be found in Figure 16. The modelled NH 3 surface concentrations are 205 compared to observation from the Dutch MAN network (Lolkema et al., 2015) and the European EMEP network (EMEP, 2016). The modelled NO 2 surface concentrations are compared to observation from AirBase (EEA, 2019).
We only used background stations. The location of these stations can be found in Figure 14.

Eddy covariance data
FLUXNET is a global network of micrometeorological towers that measure biosphere-atmosphere exchange fluxes 210 using the eddy covariance (EC) method. We used half-hourly observations from the FLUXNET2015 dataset (Pastorello et al., 2017) to validate the z 0 values. We use observations of the mean wind speed, the friction velocity, the sensible heat flux, precipitation and air temperature T a to determine z 0 values . The sites used in this study are shown in Table 1.

Updated z 0 maps
The updated z 0 maps are a composition of z 0 values derived using different methods. We distinguish three different main approaches: 1) z 0 values that depend on forest canopy height, 2) z 0 values that depend on the NDVI and 3) new z 0 values for urban areas that depend on the population density map. In addition to these three approaches, the z 0 values of some urban classes were set to new default values. An overview of the datasets that are used for each 220 DEPAC land use class is given in Table 2.
The MODIS NDVI, the MODIS LAI and the GLAS forest canopy height had to be pre-processed and homogenized in order to obtain consistent input maps that can be read into the LOTOS-EUROS model. To achieve this, we created input maps for each DEPAC class on the coordinate grid of the CORINE/Smiatek land use map in LOTOS-

225
First of all, the original datasets were re-projected to geographic coordinates. The following approach is used to deal with the different horizontal resolutions of the datasets: We used the CLC2012 map, having the highest horizontal resolution, as a basis for the computation of the updated z 0 values. For each of the other datasets, we first computed the percentages of each CORINE land cover class within every pixel. We define homogeneous pixels, consisting of nearly one CORINE land cover class, for which we will use a threshold value of 85% of the pixels area. Then we 230 isolate and use only these (nearly) homogeneous pixels to compute z 0 values for each CORINE land cover class. The methods that were applied are described in the subsequent section.

Forest canopy height derived z 0 values
The forest canopy height dataset derived from GLAS LiDAR observations is used to compute the z 0 values for each CLC2012 forest land cover class (broad-leaved forest, coniferous forest and mixed forest) that corresponds to one of 235 the DEPAC forest land use classes (4: coniferous forest and 5: deciduous forest). Several publications relate vegetation height to z 0 (e.g. (Raupach, 1994;Plate, 1971;Brutsaert, 2013)). Here we used the often used equation from (Brutsaert, 2013): The vegetation height is the most important parameter influencing turbulence near the surface, and for this reason, 240 the used parameterisation gives a reasonable estimate of z 0 , even though it ignores many other aspects that influence z 0 (e.g. density, vertical distribution of foliage). Multiple studies have illustrated that there is a seasonal variation in z 0 /h for different types of forests (e.g. (Yang and Friedl, 2003;Nakai, 2008)). The z 0 of deciduous trees is therefore additionally linked to the leaf area index to account for changes in tree foliage. The following formula is used to compute the monthly z 0 value for each deciduous forest pixel: Here the maximum roughness length z , is set to the LiDAR-derived z 0 from (eq. 4). The minimum roughness length z , represents the z 0 of leafless deciduous trees. Following the dependence of z 0 /h on LAI presented in Nakai (2008) and Yang and Friedl (2003), we set the z , to 80% of z , . Table 3 gives an overview of several studies that relate the z 0 value to the NDVI. The functions are derived for different vegetation types under specific conditions. Equations 6 to 12 are derived for different types of agricultural land. These equations are all within a reasonable range from one another for NDVI values below ~0.8. Therefore, we chose to use the average function of eq. 6 to eq. 11 to compute z 0 values for all CLC subcategories of the following DEPAC classes: "arable", "other" and "permanent crops". Figure 1 is a visualization of eq. 6 to 11 and the 255 used mean function. Figure S1 shows histogram of all NDVI values in our study area in 2014. We computed that 7.4% of all NDVI values have a NDVI > 0.8, 1.3% have a NDVI > 0.85 and only 0.04% have a NDVI > 0.9.

NDVI derived z 0 values
Virtually all NDVI values thus fall within the range where the average function does not differ much from the individual functions. The z 0 value of grasslands is in general lower than other vegetation types. The last equation, eq. 12, is specifically derived for grassland and is therefore used for all CLC subcategories that fall under the 260 DEPAC class "grass".

z 0 values for urban areas
The default z 0 of urban areas in LOTOS-EUROS was set to 2 meters. We have updated the z 0 values for urban areas to avoid possible overestimation of z 0 in sparsely populated urban areas. The updated z 0 values for CLC2012 class 1 and 2, '1: continuous urban fabric' and '2: discontinuous urban fabric' are time-invariant and coupled to the EEA 265 population density map. The z 0 values are set to 2 metres in areas with a population density higher than 5000 inhabitants/km 2 and to 1 metre in areas with a population density lower than 5000 inhabitants/km 2 . The z 0 values of the other urban subcategories, CLC2012 class 3 to 9, are updated to the proposed values for CLC classes in (Silva et al., 2007). Figure S2 shows the resulting updated z 0 values for urban areas.

270
After the computation of the z 0 values, the maps for each CORINE land cover class were merged and converted into DEPAC classes using a pre-defined conversion table. As multiple CORINE land cover may translate to one single DEPAC class, the weighted average based on the respective percentage of each CORINE land cover class was computed for each pixel. We then used linear interpolation to obtain continuous z 0 maps for each DEPAC class.
Finally, the maps were re-gridded unto the CORINE/Smiatek grid and saved into one file per month.

275
The default parameterization of the LAI in LOTOS-EUROS is replaced by the MCD15A2H LAI product from MODIS. First, we applied a coordinate transformation to obtain the data in geographical coordinates. The data was then re-gridded unto the grid of the CORINE/Smiatek land use map using linear interpolation. The quality tags were evaluated to identify pixels with default fill values from the MCD15A2H product. These fill-values were replaced by the default LAI values in LOTOS-EUROS, to avoid modelling discrepancies resulting from sudden jumps in LAI values. Finally, the values were sorted per DEPAC land use class and individual fields were created for each class as new input for LOTOS-EUROS.

z 0 values from EC measurements
We used the regression method (e.g. Graf et al., 2014. Chen et al., 2015 to compute z 0 from several eddy covariance sites. A description of the methodology and the data processing is given in this section. The wind profile in the 285 surface layer can be approximated by: here, is the measurement height, is the displacement height, is the aerodynamic roughness length, is the Von-Karman constant (=0.4), ( ) is the average wind speed, * is the friction velocity and Ψ is the integrated universal momentum function, also known as the stability correction term. Ψ is a function of , the Monin-

290
Obukhov length, which is defined as (e.g. Erisman and Duyzer, 1991): where is the air temperature, the air density (= 1.2 kg m -3 ), the heat capacity at constant pressure (=1005 J kg -1 K -1 ), the acceleration due to gravity, and the sensible heat flux. Stability correction term Ψ is in principle a non-linear function, however, for a certain stability range it can be approximated by a linear function. It is shown 295 that for moderately stable conditions (0 < < 1) stability correction term Ψ holds the following form: where is a constant. We consider a simple linear regression with offset parameter and slope parameter . If we assume that Ψ is linear, we can rewrite Eq. 1 in the following form: Now provides an estimate of ln( − )/ , and we can directly compute from ( − )/exp ( ). We use observations from 2014 only, unless stated otherwise in Table 1. For forest we assume that = (2/3) * ℎ (Maurer et al., 2013), and we use the forest canopy height derived from GLAS. For short vegetation we assume that displacement height is negligible, that is, = 0 . Graf et al., 2014 illustrated that the linearity approximation of Ψ is valid for small negative values of ( − )/ , so we first select all points where -0.1 < ( − )/ < 1. We 305 filter out observation during rainfall and where * < 0.15, as presented in Chen et al., 2015. We split our data into a group with stable conditions ( > 0) and with unstable conditions ( < 0). We assume that the is more or less constant over a period of 5 days. For each 5-day period we plot ( )/ * against ( − )/ and fit a simple line function using linear least-squares. The z 0 values are then computed from offset parameter . We compute the mean, median, standard deviation and the range of the all computed z 0 values in one year. If the computed z 0 values for stable and unstable conditions in one 5-day period differ more than 50% from their arithmetic mean they are filtered out.

Comparison of the default and updated z 0 values
We used the CORINE/Smiatek land use map to combine all the updated z 0 values into one single map. The resulting 315 composite map has a horizontal resolution of 500 by 500 metres and is shown in Figure 2.
The mean relative difference (MRD) between the default and updated z 0 values is presented in Figure 3. The largest positive differences occur in forested areas, meaning that the default z 0 values are lower than the updated z 0 values.
The largest negative deviations occur in urban areas and areas with "grass". The updated z 0 values are generally lower in the Netherlands and Belgium, and mostly higher in Germany. Table 4 gives an overview of the default z 0 320 values in LOTOS-EUROS and the mean and standard deviation of the new z 0 values for each of the DEPAC land use classes. The updated z 0 values for "arable land", "coniferous forest", "deciduous forest" and "other" are on average higher than the default z 0 values in LOTOS-EUROS. The updated z 0 values for "grass", "permanent crops" and "urban" are on average lower than the default z 0 values in LOTOS-EUROS.

325
We compared the updated z 0 values to z 0 values from several studies (Wieringa, 1993;Silva et al., 2007;Troen and Petersen, 1989;Lankreijer et al., 1993;Yang and Friedl, 2003), and z 0 values used in other CTMs (Simpson et al., 2012;Mailler et al., 2017). Table 5 gives an overview of these z 0 values. There is in general good agreement with these z 0 values, and the updated z 0 values mostly fall within the stated ranges. The updated z 0 values for coniferous and deciduous forest are on the high side compared to these studies. A histogram of the forest canopy heights 330 derived from GLAS within our study area is given in Figure S1. These differences can in part be explained by the occurrence of relatively tall forest canopy (~30 meters) in the dataset, especially in forest in southern Germany, whereas most of these studies either assumed or studied shorter trees. Another possible explanation lies in the fact that we used a relatively large conversion factor of 0.136 (eq. 4), whereas a factor of 0.10 is also used quite often.

365
deciduous broadleaf forest sites reach the highest LAI values in the growing season. There is less variation in the LAI for evergreen needle leaf forest sites. However, the wintertime LAI values seem to be unrealistically low.

Implications for modelled N r deposition fields
In the following section, the impact of the updated LAI and z 0 values on modelled N r deposition fields in LOTOS-EUROS is discussed. A total of four different LOTOS-EUROS runs are compared to examine the individual effect 370 of the updated LAI and z 0 values on the modelled N r distributions and deposition fields. The first run, named "LE default ", is the standard run using default LAI and z 0 values. The second run, named "LE LAI ", uses updated LAI values and the default z 0 values. The third run, named "LE z0 ", uses updated z 0 values and the default LAI values.
The fourth run, named "LE z0+LAI ", uses both updated LAI and z 0 values. From now on, we will refer to the outputs of these different runs using the abovementioned abbreviations. Figure S3 shows the division of the total terrestrial N r deposition over Germany, the Netherlands and Belgium into different N r compounds for each of the model runs. A relatively larger portion of the total N r deposition is attributed to oxidized forms of N r in Germany. The reduced forms of N r predominate in the Netherlands and Belgium. The largest change in total N r deposition occurs in Belgium (+6.19%) due to the inclusion of the MODIS LAI. This 380 corresponds to the relative increase in LAI values here. The inclusion of the updated z 0 values lead to a minor decrease in total N r deposition in the Netherlands (-1.45%) and Belgium (-1.13%), and a minor increase in Germany (+0.44%).

Effect on wet and dry N r deposition
We examined the direct effect of the updated LAI and z 0 values on the modelled dry N r deposition, as well as the 385 related indirect effect in modelled wet N r deposition. Figure 10 shows the dry and wet N r deposition in kg N ha -1 in 2014, modelled with the updated LAI and z 0 values as input in LOTOS-EUROS. Figure 11 shows the relative changes in the total amount of dry and wet N r deposition of the different runs with respect to the default run. The combined effect shows an increase of the amount of dry N r deposition over most parts of Belgium and Germany.
The amount of wet N r deposition decreases over most parts of Germany and eastern Belgium, but remains 390 unchanged in north-western parts of Germany. We observe a decrease in total N r deposition in the Netherlands. In general, we observe changes ranging from approximately -20% to +30% in the total amount of dry N r deposition.
The changes in wet N r deposition are smaller in magnitude and range from approximately -3% to +3%.

Effect on reduced and oxidized N r deposition
We split up the total Nr deposition into NH x (NH 3 and NH 4 + ) and NO y (NO and NO 2 and NO 3 and HNO 3 ) deposition, to look at the effect of the updated LAI and z 0 input maps on the deposition of reduced and oxidised forms of N r , respectively. Figure 12 shows the modelled NH x and NO y deposition in kg N ha -1 in 2014, including the updated LAI and z 0 input values. Figure 13 shows the relative changes (%) in the total NH x and NO y of the different runs with respect to the default run of LOTOS-EUROS. The updated z 0 values have a larger impact on the NH x deposition than on the NO y deposition. The implementation of the updated z 0 values has led to a decrease in NH x 400 deposition over most of the Netherlands, and western Belgium, driven by the large fraction of grassland here. The updated LAI values led to relatively more NH x deposition in Belgium. The updated LAI values led to an increase of NO y deposition in almost all areas within the modelled region, except for some urban areas. Moreover, the impact seems to be limited in large forests located in background areas. Table 7 gives an overview of changes in the distribution of the land use specific fluxes in the whole study area (Germany, the Netherlands and Belgium combined) for the different runs. The most distinct changes in N r deposition are due to the updated LAI values ("LE LAI "), where we observe an increase in total N r deposition on urban areas (+ 16.62%) and arable land (+ 9.53%), and a decrease on coniferous forests (-9.36%). This coincides with the categories where we observe the largest changes in LAI values. The default LAI values in urban areas were first set to zero for all urban DEPAC categories. The MODIS LAI values, however, are only zero in densely populated areas and areas with industry. The main effects of the updated z 0 values ("LE z0 ") can be observed for grass (-3.95 %), permanent crops (+ 3.27) and arable land (-3.17 %). In the combined run, "LE z0+LAI ", we observe an amplified effect in total N r deposition over grass (-8.05%) and arable land (+ 12.88%). The impact of the individual parameters on the remaining land use categories is attenuated in this run.

Implications for N r distributions
The changes in N r deposition amounts induce an effect in the modelled distribution of nitrogen components. Here, we look at the effect of the updated LAI and z 0 values on NH 3 and NO 2 surface concentrations. Figure 14 shows

Comparison to in-situ measurements
To analyse the effect of the updated LAI and z 0 values, we compared our model output to the available in-situ observations. Due to the lack of available dry deposition measurements, we decided to use NH 4 + and NO 3 wet deposition and NH 3 and NO 2 surface concentrations measurements instead. The distribution of the wet deposition stations is shown in Figure 16, as well as the modelled mean NH 4 + (left) and NO 3 -(right) wet deposition in 2014.

435
The locations of the stations that measure the NH 3 and NO 2 surface concentrations are shown in Figure 14.
The relationships between the modelled and observed fields are evaluated using the Pearson's correlation coefficient (r), the RMSD and the coefficients (slope, intercept) of simple linear regression. Table S1 shows these measures for the comparison with monthly mean NO 3 wet deposition, NH 4 + wet deposition, and the monthly mean NH 3 and NO 2 surface concentrations in 2014. Table S1 shows the same statistics but computed per DEPAC land use class.

440
Overall, we do not observe large changes in the shown measures due to the inclusion of the updated LAI and z 0 values on a yearly basis. The model underestimates NO 3 wet deposition, and NH 4 + to a lesser extent, too. The modelled NH 3 surface concentrations are in general higher than observed concentrations. The opposite is true for NO 2 , here, the modelled surface concentrations are lower than the observed concentrations. The computed measures did not change drastically due to the inclusion of the updated z 0 and LAI values. Figure S4 shows the monthly mean NO 3 wet deposition, NH 4 + wet deposition, NH 3 surface concentration and NO 2 surface concentrations for the different model runs and the mean of the corresponding in situ observations. The Both Table S1, Table S2 and Figure 16 illustrate that the comparability of the modelled wet deposition and surface concentration fields to the available in-situ measurements did not change significantly. The impact of the updated LAI and z 0 values on these fields is largely an indirect effect of the more distinct changes in the dry deposition, and thus too small to lead to any drastic changes. We conclude that we are unable to demonstrate any major 460 improvements with the use of the currently available in-situ measurements.

Discussion
This paper aimed to improve the surface characterization of LOTOS-EUROS through the inclusion of satellite- The z 0 is closely related to the geometric features and distributions of the roughness elements in a certain area. The 475 updated z 0 values are linked to specific land use pixels and are therefore assumed reasonable estimates for moderately homogeneous areas with this specific land use type. There are various approaches to combine these z 0 values into an 'effective' roughness for larger, mixed areas (e.g. (Claussen, 1990;Mason, 1988)). The LOTOS-EUROS model uses logarithmic averaging to compute an effective roughness for an entire model pixel. This averaging step seems to be one of the reasons why the effect of our updated z 0 values on the deposition fields is 480 limited. To illustrate this, the relative change in total dry NH 3 deposition due to the updated z 0 values were computed and shown in Figure 17. We used increasing threshold percentages to sort the NH 3 deposition on a model pixel level per land use type and fraction. Figure 17 shows that the differences in total NH 3 deposition between the two runs increase with increasing land use fraction. The model pixels that mostly consist of one type of land use seem to show the largest change in NH 3 deposition. The change thus appears to be less distinct in pixels that have a 485 higher degree of mixing. Most of the model pixels largely contain mixtures of different land uses on the current model resolution. As a result, averaging of z 0 on a model pixel level is thus likely to cause a levelling effect on the current model resolution. The impact of the updated z 0 values is therefore expected to be larger at a higher model resolution. The use of another approach for computing the 'effective' roughness could potentially lead to stronger changes in the modelled deposition fields.

490
Moreover, we should also consider the limitations of the datasets used in this study. The forest canopy height map used in this study has been validated against 66 FLUXNET sites (Simard et al., 2011). The results showed that RMSE = 4.4 m and R 2 = 0.7 after removal of 7 outliers. For the FLUXNET forest sites used in this study, we compared the forest canopy heights from GLAS to the maximum forest canopy height at the FLUXNET sites taken from (Flechard et al., 2019). For all but one site (DE-Hai), the forest canopy heights from GLAS were lower than 495 this value (Table S3). This method could potentially be improved by using another product with either a higher precision or resolution. For modelling studies on a national level one could or instance consider the use of airborne LiDAR point clouds to retrieve forest canopy heights. This procedure, although it is computationally expensive, would allow us to create high resolution z 0 maps.
The MODIS-LAI at the FLUXNET sites showed realistic seasonal variations for most land use classifications, 500 except for relatively low winter-time values for evergreen needle leaf forests. The previous versions of the MODIS-LAI have been validated in many studies (e.g. (Fang et al., 2012;Wang et al., 2004;Kobayashi et al., 2010)), showing an overall good agreement with ground observed LAI values and other LAI products. The seasonality in LAI is properly captured for most biomes, but unrealistic temporal variability is observed for forest due to infrequent observations. Also, the previous versions overestimate LAI for forests (Fang et al., 2012;Kobayashi et al.,  Some studies (e.g. Tian et al., 2004) have reported an underestimation of the MODIS-LAI in presence of snow-510 cover, particularly affecting evergreen forests. There was only a limited amount of snowfall in our study region, so this did not lead to any problems. However, this issue should be carefully considered when using the MODIS-LAI for regions with frequent snow-cover, like Scandinavia. LOTOS-EUROS uses meteorological data to determine what regions are covered with snow. If a certain location is snow-covered, the standard parameterization for the canopy resistance is not used. LOTOS-EUROS uses a pre-defined value for the canopy resistance instead. As such, these low MODIS-LAI values do not affect the modelled deposition in LOTOS-EUROS during snow-cover.
Though the issues with the MODIS-LAI should be considered with care, the spatial and temporal distribution of these LAI values is more realistic than that of the default LAI values used in LOTOS-EUROS. The same holds for the updated, time-variant z 0 values. The representation of the growing season is now more realistic due to their dependence on NDVI and LAI values. Moreover, the z 0 values for forests now also have a clear spatial variation, 520 such as a latitudinal gradient with increasing z 0 values towards the south of Germany. These type of patterns can simply not be captured by fixed values.
We evaluated the effect of updated z 0 and LAI values on modelled N r distribution and deposition fields. The distribution of the relative changes in deposition of the reduced and oxidized forms of reactive nitrogen showed a similar pattern. Here, the updated z 0 values led to a variation of ~±8%, and the updated LAI values led to variations 525 of ~±30% in both fields. The dry deposition fields were most sensitive to changes in z 0 and LAI, as these varied from approximately -20% to +20% with the updated z 0 map, and from -20% even up to +30% with the MODIS LAI values. As a result, we observed a shift from wet to dry deposition, except for the Netherlands, where we observe an opposite shift, from dry to wet deposition. Moreover, we observed a redistribution of N r deposition over different land use classes on a sub grid level. To illustrate the potential consequences on a local scale, we computed the 530 critical load exceedances for deciduous and coniferous forest ( Figure 18) using critical loads of 10 kg following Bobbink et al. (2010b). Compared to the default run, the changes may be sizable locally, ranging from approximately -3 kg up to +2 kg for deciduous forest and even over -3 kg for coniferous forest.
The uncertainties of the LAI and z 0 input data are but one aspect of the model uncertainty of CTMs. The model uncertainty has several other origins, like the physical parameterizations (e.g. deposition velocities) and the 535 numerical approximations (e.g. grid size). Two of the most important uncertainties related to deposition modelling are the emissions and the surface exchange parameterization. The emission inventories for reactive nitrogen hold a relatively large uncertainty. The uncertainty of the reported annual total NH 3 emissions is estimated to be at least ±30%. This is mainly due to the diverse nature of agricultural emission sources, leading to large spatio-temporal variations. The annual NO x emissions total hold a lower uncertainty, of approximately ±20% (Kuenen et al., 2014).

540
Emissions at specific locations, especially for NH 3 , are even more uncertain due to assumptions made in the redistribution and timing of emissions. A recent paper of Dammers et al., 2019, for instance, found that satellitederived NH 3 emissions of large point sources are a factor 2.5 higher than those given in emission N-inventories.
The surface-exchange parameterization is another source of uncertainty. The complexity of the NH 3 surface exchange schemes in CTMs is usually low compared to the current level of process understanding (Flechard et al., 545 2013). Moreover, large discrepancies exist between deposition schemes. Flechard et al. (2011) for instance showed that the differences between four dry deposition schemes for reactive nitrogen can be as large as a factor 2-3.
This work has shown that changes in two of the main deposition parameters (LAI, z 0 ) can already lead to distinct, systematic changes (~30%) in the modelled deposition fields. This demonstrates the model's sensitivity toward these input values, especially the LAI. In addition to the known uncertainty involved with surface exchange parameterization itself, this further stresses the need for further research. Another important aspect that should receive more attention is the validation of the dry deposition fields with in-situ dry deposition measurements. Here we illustrated the need for direct validation methods, as relatively large changes in modelled dry deposition field cannot be verified sufficiently by surface concentration and wet deposition measurements.
The surface-atmosphere exchange remains one of the most important uncertainties in deposition modelling. The use 555 of satellite products to derive LAI and z 0 values can help us to represent the surface characterization in models more accurately, which in turn might help us to minimize the uncertainty in deposition modelling. The approach to derive high resolution, dynamic z 0 estimates presented here can be linked to any land use map and is as such transferable to many different models and geographical areas.         (7) 16.62 -0.37 6.53 Other (6,8,9) 3.45 0.58 3.31