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 (N_r) 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 N_r flows. In these models, the parameterization of the atmosphere–biosphere exchange of N_r 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 (z₀) input maps. As LAI, we use the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD15A2H product. The monthly z₀ input maps presented in this paper are a function of satellite-derived 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 z₀ values are updated, too, and linked to a population density map. The approach to derive these dynamic z₀ estimates can be linked to any land use map and is as such transferable to other models. We evaluated the sensitivity of the modelled N_r deposition fields in LOng Term Ozone Simulation – EURopean Operational Smog (LOTOS-EUROS) v2.0 to the abovementioned changes in LAI and z₀ inputs, focusing on Germany, the Netherlands and Belgium. We computed z₀ values from FLUXNET sites and compared these to the default and updated z₀ 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 z₀) to 0.60 (updated z₀). The implementation of these updated LAI and z₀ input maps led to local changes in the total N_r 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.

1 Introduction

The nitrogen (N) cycle has been continuously disrupted by human activity over the past century (Fowler et al., 2015; Galloway et al., 2004, 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₂, have manifested in the environment contributing to acidification and eutrophication of sensitive terrestrial and aquatic ecosystems (Bobbink et al., 2010a; Paerl et al., 2014). NO_x and NH₃ 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 (Erisman et al., 2011; 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 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 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₃ fields are in general uncertain due to the highly reactive nature and the uncertain lifetime of NH₃ in the atmosphere. More importantly, NH₃ emissions that are used as model input are very complex to estimate and remain highly uncertain (Reis et al., 2009; Behera et al., 2013), for example, due to the diversity in NH₃ volatilization rates originating from different agricultural practices. Recently, a lot of effort has been made to improve the spatiotemporal distributions of bottom-up NH₃ emissions (e.g. Hendriks et al., 2016; 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 has already enabled us to create large-scale emission estimates for various pollutants (e.g. Curier et al., 2014; Abida et al., 2017), for instance, for NO₂ and will likely also be used for large-scale NH₃ 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 often pragmatic and computed with relatively simple empirical functions (e.g. following Wesely, 1989; Emberson et al., 2000; Erisman et al., 1994), mostly linked to land use classification maps. The parameterization of the atmosphere–biosphere exchange of N_r components that is used in models is largely based on surface exchange measurements and is therefore very uncertain, especially for NH₃ (Schrader et al., 2016). The deposition strengths in models may vary tremendously depending on the used deposition parameterization and velocities (Wu et al., 2018; Schrader and Brummer, 2014; Flechard et al., 2013). Moreover, inter-model discrepancies in deposition fluxes may also arise from differences in the used input variables. Here, we focus on the leaf area index (LAI) and the roughness length (z₀) input values. The deposition velocity is often parameterized using both the LAI and the z₀. Currently, most models use fixed, land-use-specific values for both parameters. In practice, however, spatial as well as seasonal variation is observed. In this paper, we aim to improve the deposition flux modelling by using more realistic, space- and time-variant LAI and z₀ 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 parameterized using the LAI. The LAI can be determined in the field using either direct 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-known example of the latter is the LAI product from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument, which we will use in this study.

The z₀ 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:

$$\begin{array}{}\text{(1)}& {\displaystyle }U\left(z\right)={\displaystyle \frac{u_{*}}{k}}\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right).\end{array}$$

U(z) represents the mean wind speed, u_* the friction velocity and k the von Kármán constant. Here, z₀ is a constant that represents the height at which the wind speed theoretically becomes zero. The z₀ can be estimated from in situ wind speed measurements using bulk transfer equations. More recently, several studies have shown that z₀ for specific, uniform land cover types can also be estimated from optical remote sensing measurements, for instance, using vegetation indices (e.g. Xing et al., 2017; Yu et al., 2016; Bolle and Streckenbach, 1993; Hatfield, 1988; Moran, 1990). The z₀ 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 large potential for improvements to the representativeness of the surface characterization in chemistry transport models. Here, we present an approach to derive monthly z₀ input maps using satellite-derived normalized difference vegetation index (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₀ values by comparing them to z₀ values computed from FLUXNET observations. We also update the z₀ values for urban areas, using a population density map. We use the updated z₀ values, as well as the MODIS LAI, as input in LOng Term Ozone Simulation – EURopean Operational Smog (LOTOS-EUROS) to illustrate the effect on transport and deposition modelling of N_r components. We evaluate the sensitivity of the N_r deposition 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 ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and surface concentration measurements of NH₃ and NO₂.

2 Model and datasets

2.1 LOTOS-EUROS

2.1.1 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 km 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 m. The next layer is a mixing layer, 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 European Centre for Medium-Range Weather Forecasts (ECMWF) meteorology. The gas-phase chemistry uses the Netherlands Organisation for Applied Scientific Research (TNO) carbon bond mechanism (CBM) IV scheme (Schaap et al., 2009) and the anthropogenic emissions from the TNO Monitoring Atmospheric Composition and Climate (MACC) III emission database (Kuenen et al., 2014). The wet deposition parameterization is based on the Comprehensive Air Quality Model with Extensions (CAMx) approach and includes both in-cloud and below-cloud scavenging (Banzhaf et al., 2012). LOTOS-EUROS makes use of the Coordination of Information on the Environment (CORINE)/Smiatek land use map to determine input values for surface variables.

2.1.2 Dry deposition

The dry deposition flux of gases is computed following the resistance approach, in which the exchange velocity V_d is equal to the reciprocal sum of the aerodynamic resistance R_a, the quasi-laminar boundary layer resistance R_b and the canopy resistance R_c:

$$\begin{array}{}\text{(2)}& V_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{b}}+R_{\mathrm{c}}}}.\end{array}$$

R_a and R_b are both influenced by the wind profile, which is computed with (Eq. 1). The wind profile, in turn, depends on roughness length z₀ associated with different land use classes. The aerodynamic resistance R_a is computed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{R}_{\mathrm{a}}=\underset{z_{\mathrm{0}}}{\overset{h}{\int }}{\displaystyle \frac{\mathrm{\Phi }\left(\frac{z}{L}\right)}{\mathit{\kappa }{u}_{*}z}}\mathrm{d}z.\end{array}$$

Here, h is the canopy height, which is pre-defined per land use class. The empirical function $\mathrm{\Phi }\left(\frac{z}{L}\right)$ is taken from Businger et al. (1971) and depends on the state of the atmosphere. Depending on the value of Monin–Obukhov length L, the following equations are used:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}& \text{For a stable atmosphere,}(L>\mathrm{0}):\\ & {\mathrm{\Phi }}_{\mathrm{s}}\left({\displaystyle \frac{z}{L}}\right)=\mathrm{1}+\mathrm{4.7}\left({\displaystyle \frac{z}{L}}\right)\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}& \text{For an unstable atmosphere,}(L<\mathrm{0}):\\ & {\mathrm{\Phi }}_{\mathrm{u}}\left({\displaystyle \frac{z}{L}}\right)={\left(\mathrm{1}-\mathrm{15}\left({\displaystyle \frac{z}{L}}\right)\right)}^{-\mathrm{0.25}}.\end{array}\end{array}$$

For a neutral atmosphere, it is equal to unity. The Monin–Obukhov length L is dependent on z₀ and is determined as follows (Majumdar and Ricklin, 2010):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{L}}=S\left(a_{\mathrm{1}}+a_{\mathrm{2}}{S}^{\mathrm{2}}\right)z_{\mathrm{0}}^{(b_{\mathrm{1}}+b_{\mathrm{2}}\left|S\right|+b_{\mathrm{3}}{S}^{\mathrm{2}})}.\end{array}$$

Here, a₁, a₂, b₁, b₂ and b₃ are constants (0.004349, 0.003724, −0.5034, 0.2310 and −0.0325, respectively) and $S=-\mathrm{0.5}$ (3.0–0.5 u_s|CE|) with near-surface wind speed u_s and exposure factor CE, depending on cloud cover and solar zenith angle. The wind speed at a reference height (10 m) is used to obtain the friction velocity u_*:

$$\begin{array}{}\text{(7)}& {\displaystyle }{u}_{*}=\mathit{\kappa }{u}_{\mathrm{10}\mathrm{m}}/\ln \left({\displaystyle \frac{z_{\mathrm{10}\mathrm{m}}}{z_{\mathrm{0}}}}\right).\end{array}$$

The quasi-laminar boundary layer resistance R_b is a function of the cross-wind leaf dimension L_d and the wind speed at canopy top, u(h), following the parameterization presented in McNaughton and Van Den Hurk (1995):

$$\begin{array}{}\text{(8)}& {\displaystyle }{R}_{\mathrm{b}}=\mathrm{1.3}\cdot \mathrm{150}\cdot \sqrt{{\displaystyle \frac{L_{\mathrm{d}}}{u\left(h\right)}}}.\end{array}$$

L_d 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, L_d, and subsequently R_b, is equal to zero. The canopy resistance R_c is computed using the Deposition of Acidifying Compounds v3.11 (DEPAC3.11) module (van Zanten et al., 2010). R_c is a parallel system of the resistances of three different pathways, the external leaf surface or cuticular resistance R_w, the effective soil resistance R_soil,eff and the stomatal resistance R_s, and is defined as

$$\begin{array}{}\text{(9)}& {\displaystyle }{R}_{\mathrm{c}}={\left({\displaystyle \frac{\mathrm{1}}{R_w}}+{\displaystyle \frac{\mathrm{1}}{R_{\text{soil,eff}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{s}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The external leaf surface resistance R_w 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 R_soil,eff is the sum of the in-canopy resistance R_inc and the soil resistance R_soil. Soil resistance R_soil has a fixed value, depending on the land use class and conditions (frozen, wet or dry). In the case of $u_{*}>\mathrm{0}$, the in-canopy resistance for arable land, permanent crops and forest is computed as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }{R}_{\mathrm{inc}}={\displaystyle \frac{\mathrm{14}h\text{SAI}}{u_{*}}}.\end{array}$$

For $u_{*}<\mathrm{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 by correction factors for photoactive radiation, temperature and vapour pressure deficit to obtain the stomatal conductance G_s. The R_s is then equal to 1∕G_s. The resistance parameterizations differ with land use type. A total of nine different land use classes are used in DEPAC. LOTOS-EUROS uses a fixed z₀ 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₃ is included in the implementation of the DEPAC3.11 module (Wichink Kruit et al., 2012), allowing emissions of NH₃ under certain atmospheric conditions. More information about the most recent version of the model can be found in Manders et al. (2017).

2.2 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₀ maps are described. Finally, the in situ observations used for validation are discussed in the last paragraphs.

2.2.1 MCD15A2H leaf area index

The satellite-derived LAI is a combined product of the MODIS instruments aboard the Terra and Aqua satellites (Myneni et al., 2015). The LAI algorithm compares bidirectional spectral reflectances observed 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 sixth version of the product, MCD15A2H, which has a temporal resolution of 8 d and a spatial resolution of 500 m.

2.2.2 MYD13A3 NDVI

NDVI is a vegetation index computed with reflectances observed by the MODIS instrument aboard the Aqua satellite (Didan, 2015). The NDVI is the normalized transform of the near infrared to the red reflectance and is expressed by

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathrm{NDVI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{red}}}{{\mathit{\rho }}_{\mathrm{NIR}}+{\mathit{\rho }}_{\mathrm{red}}}}.\end{array}$$

We used the MYD13A3 product, which is the monthly NDVI product with a spatial resolution of 1 km.

Table 1 Overview of eddy covariance sites used to compute z₀. The following abbreviations for land use types are used: DBF – deciduous broadleaf forest, ENF – evergreen needleleaf coniferous forest, MF – mixed forest, CRO – croplands, WET – permanent wetlands, GRA – grasslands.

References for each Site ID: BE-Bra: Janssens (2016); BE-Lon: De Ligne et al. (2016a); BE-Vie: De Ligne et al. (2016b); DE-Akm: Bernhofer et al. (2016a); DE-Geb: Brümmer et al. (2016); DE-Gri: Bernhofer et al. (2016b); DE-Hai: Knohl et al. (2016); DE-Kli: Bernhofer et al. (2016c); DE-Obe: Bernhofer et al. (2016d); DE-RuR: Schmidt and Graf (2016a); DE-RuS: Schmidt and Graf (2016b); DE-Seh: Schneider and Schmidt (2016); DE-SfN: Klatt et al. (2016); DE-Tha: Bernhofer et al. (2016e); NL-Hor: Dolman et al. (2016); NL-Loo: Moors and Elbers (2016); FR-Fon: Berveiller et al. (2016); CH-Cha: Hörtnagl et al. (2016a); CH-Fru: Hörtnagl et al. (2016b); CH-Lae: Hörtnagl et al. (2016c); CH-Oe1: Ammann (2016); CH-Oe2: Hörtnagl et al. (2016d); CZ-wet: Dušek et al. (2016).

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

2.2.4 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 spatial resolution of 100 m.

2.2.5 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 latest version of the product, CLC2012, in this study.

Table 2 An overview of the datasets that are used to derive z₀ input values for each DEPAC land use category.

2.2.6 In situ reactive nitrogen observations

The modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ wet deposition fluxes are compared to observations of wet-only samplers. We used observation from the Dutch Landelijk Meetnet Regenwatersamenstelling (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 Fig. 16. The modelled NH₃ surface concentrations are compared to observation from the Dutch MAN network (Lolkema et al., 2015) and the European Monitoring and Evaluation Programme (EMEP) network (EMEP, 2016). The modelled NO₂ surface concentrations are compared to observation from AirBase (EEA, 2019). We only used background stations. The location of these stations can be found in Fig. 14.

2.2.7 Eddy covariance data

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

Table 3 Studies that relate the aerodynamic roughness length (z₀) to satellite-derived NDVI values for specific land cover types and conditions.

Figure 1 Several functions that relate the roughness length z₀ (m) to the normalized difference vegetation index (NDVI). The dotted line shows the average of all the functions. This function is used to compute NDVI-dependent z₀ values for the subcategories within DEPAC classes “arable”, “other” and “permanent crops”.

3 Methodology

3.1 Updated z₀ maps

The updated z₀ maps are a composition of z₀ values derived using different methods. We distinguish three different main approaches: (1) z₀ values that depend on forest canopy height, (2) z₀ values that depend on the NDVI and (3) new z₀ values for urban areas that depend on the population density map. In addition to these three approaches, the z₀ values of some urban classes were set to new default values. An overview of the datasets that are used for each 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-EUROS.

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₀ 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 isolate and use only these (nearly) homogeneous pixels to compute z₀ values for each CORINE land cover class. The methods that were applied are described in the subsequent section.

3.1.1 Forest-canopy-height-derived z₀ values

The forest canopy height dataset derived from GLAS lidar observations is used to compute the z₀ values for each CLC2012 forest land cover class (broadleaved forest, coniferous forest and mixed forest) that corresponds to one of the DEPAC forest land use classes (4: coniferous forest and 5: deciduous forest). Several publications relate vegetation height to z₀ (e.g. Raupach, 1994; Plate, 1971; Brutsaert, 2013). Here, we used the often-used equation from Brutsaert (2013):

$$\begin{array}{}\text{(12)}& z_{\mathrm{0}}=\mathrm{0.136}\cdot h.\end{array}$$

The vegetation height is the most important parameter influencing turbulence near the surface, and for this reason, the used parameterization gives a reasonable estimate of z₀, even though it ignores many other aspects that influence z₀ (e.g. density, vertical distribution of foliage). Multiple studies have illustrated that there is a seasonal variation in z₀∕h for different types of forests (e.g. Yang and Friedl, 2003; Nakai, 2008). The z₀ 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₀ value for each deciduous forest pixel:

$$\begin{array}{}\text{(13)}& z_{\mathrm{0}}\left(\mathrm{LAI}\right)=z_{\mathrm{0},\min }+{\displaystyle \frac{\mathrm{LAI}-\mathrm{LAI},\min }{\mathrm{LAI},max-\mathrm{LAI},\min }}(z_{\mathrm{0},\max }-z_{\mathrm{0},\min }).\end{array}$$

Here, the maximum roughness length z_0,max is set to the lidar-derived z₀ from (Eq. 12). The minimum roughness length z_0,min represents the z₀ of leafless deciduous trees. Following the dependence of z₀∕h on LAI presented in Nakai (2008) and Yang and Friedl (2003), we set the z_0,min to 80 % of z_0,max.

Figure 2 Composite map of the new z₀ values. The yearly mean is displayed for land use classes with time-variant z₀ values.

Figure 3 Mean relative difference (%) of the updated z₀ values with respect to the default z₀ values.

Table 4 An overview of the default and the adjusted roughness length of each DEPAC class used in LOTOS-EUROS for Germany, the Netherlands and Belgium. The datasets that are used to derive the updated z₀ values are given in the last column. The mean z₀ value of the DEPAC classes computed using the MODIS NDVI product is the yearly mean of all monthly z₀ values.

Table 5 Roughness length values from different types of studies. The first column states the global land use category of the z₀ values. The second column states the (range of) z₀ values, as well as the specific type of land use they are derived for. The third column shows the reference. ¹ Literature study. ² Model input (EMEP Meteorological Synthesizing Centre – West; MSC-W). ³ Model input (CHIMERE).

3.1.2 NDVI-derived z₀ values

Table 3 gives an overview of several studies that relate the z₀ value to the NDVI. The functions are derived for different vegetation types under specific conditions. Equations (12)–(20) 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 Eqs. (12)–(19) to compute z₀ values for all CLC subcategories of the following DEPAC classes: “arable”, “other” and “permanent crops”. Figure 1 is a visualization of Eqs. (12)–(19) and the used mean function. Figure S1 in the Supplement shows a 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. Virtually all NDVI values thus fall within the range where the average function does not differ much from the individual functions. The z₀ value of grasslands is in general lower than other vegetation types. The last equation, Eq. (20), is specifically derived for grassland and is therefore used for all CLC subcategories that fall under the DEPAC class “grass”.

3.1.3 z₀ values for urban areas

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

3.2 LAI and z₀ in LOTOS-EUROS

After the computation of the z₀ values, the maps for each CLC class were merged and converted into DEPAC classes using a pre-defined conversion table. As multiple CORINE land covers may translate to a 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₀ maps for each DEPAC class. Finally, the maps were regridded unto the CORINE/Smiatek grid and saved into one file per month.

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 were then regridded 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.

4 Results

4.1 Comparison of the default and updated z₀ values

We used the CORINE/Smiatek land use map to combine all the updated z₀ values into a single map. The resulting composite map has a horizontal resolution of 500 by 500 m and is shown in Fig. 2.

The mean relative difference (MRD) between the default and updated z₀ values is presented in Fig. 3. The largest positive differences occur in forested areas, meaning that the default z₀ values are lower than the updated z₀ values. The largest negative deviations occur in urban areas and areas with “grass”. The updated z₀ values are generally lower in the Netherlands and Belgium, and mostly higher in Germany. Table 4 gives an overview of the default z₀ values in LOTOS-EUROS and the mean and standard deviation of the new z₀ values for each of the DEPAC land use classes. The updated z₀ values for “arable land”, “coniferous forest”, “deciduous forest” and “other” are on average higher than the default z₀ values in LOTOS-EUROS. The updated z₀ values for “grass”, “permanent crops” and “urban” are on average lower than the default z₀ values in LOTOS-EUROS.

Table 6 Comparison of the computed z₀ values from FLUXNET observations and the corresponding satellite-derived z₀ values. The table shows the mean and standard deviation of all z₀ values in 1 year. For forest, only the maximum z₀ value is given.

Figure 4 Comparison of the updated z₀ values (x axis) to the z₀ values derived from EC measurements (y axis). The error bars indicate the standard deviation.

4.2 Comparison to z₀ values from other studies

We compared the updated z₀ values to z₀ values from several studies (Wieringa, 1993; Silva et al., 2007; Troen and Petersen, 1989; Lankreijer et al., 1993; Yang and Friedl, 2003) and z₀ values used in other CTMs (Simpson et al., 2012; Mailler et al., 2017). Table 5 gives an overview of these z₀ values. There is in general good agreement with these z₀ values, and the updated z₀ values mostly fall within the stated ranges. The updated z₀ values for coniferous and deciduous forest are on the high side compared to these studies. A histogram of the forest canopy heights derived from GLAS within our study area is given in Fig. S1 in the Supplement. These differences can in part be explained by the occurrence of relatively tall forest canopy (∼30 m) in the dataset, especially in forests 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. 12), whereas a factor of 0.10 is also used quite often.

Figure 5 Seasonal variation of the computed z₀ values at the FLUXNET cropland sites in 2014 (left axis) and the corresponding z₀ values estimated from NDVI values (right axis). The green lines indicate the estimated z₀ values derived from FLUXNET measurements and the purple lines indicate the z₀ values derived from NDVI values. The assumptions used to compute the z₀ values are shown in the titles of the figures.

Figure 6 The yearly mean MODIS leaf area index in 2014.

Figure 7 The absolute differences (LAI_MODIS – LAI_default) between the MODIS leaf area index and the default leaf area index in LOTOS-EUROS.

Figure 8 Seasonal variation of the default and MODIS LAI values per DEPAC class. The black line represents the mean MODIS LAI of all pixels within the modelled grid for that particular DEPAC class; the ranges represent the mean plus and minus the standard deviation of the MODIS LAI. The red line depicts the default LAI values in LOTOS-EUROS.

4.3 Comparison to z₀ values derived from EC measurement sites

We computed the z₀ values of the EC sites, and we refer to Sect. S1 in the Supplement for a description of the methodology that is used to derive these values. We compared the z₀ values based on their land use stated by FLUXNET to avoid issues arising from discrepancies in land use classifications. The forest sites (DBF, ENF and MF) are compared to z₀ values derived from GLAS. The cropland and wetland sites (CRO and WET) are compared to the NDVI-dependent z₀ values derived using the mean function shown in Fig. 1. The grassland sites (GRA) are compared to the NDVI-dependent z₀ values for grassland specifically. The results per site are given in Table 6. Figure 4 shows the comparison of the z₀ values from EC measurements and the updated z₀ values for different land use classes. The z₀ values for forests match quite well. Differences between the two can in part be explained by the influence of topography, which is not accounted for in the z₀ values that are derived using forest canopy height only. The z₀ values for short vegetation seem to be overestimated for crops and underestimated for grassland and wetland sites. The underestimation of some grassland and wetland sites can be explained by the large inter-site differences in vegetation cover. Some of the FLUXNET sites classified as grasslands are, for instance, mostly covered with short grass only (for instance, Oensingen), while other sites are covered with relatively tall herbaceous vegetation, such as reeds (for instance, Horstermeer). Compared to the default z₀ values in LOTOS-EUROS, the root mean square difference (RMSD) improved from 0.76 (default z₀) to 0.60 (updated z₀). For forest, the RMSD improved from 1.23 (default z₀) to 0.96 (updated z₀). For short vegetation, the RMSD also decreased, from 0.22 (default z₀) to 0.19 (updated z₀). Figure 5 shows the comparison of the seasonal variation in computed and satellite-derived z₀ values for the FLUXNET sites classified as crops in 2014. We can once more observe a clear offset between the two. The FLUXNET z₀ values go to near-zero values in wintertime, whereas the satellite-derived z₀ values never drop below 0.12 m. This seems to be due to the distribution of the NDVI values (Fig. S1), which shows that the NDVI >0.4 most of the time. The seasonal patterns, on the other hand, seem to match well in general, even though the satellite-derived z₀ values rise somewhat earlier in the year.

Figure 9 Seasonal variation of the MODIS LAI at individual FLUXNET sites with different land use classifications.

Figure 10 The modelled amount of dry (a) and wet (b) deposition in kg N ha⁻¹ in 2014.

4.4 Comparison of the default and MODIS LAI

The yearly mean MODIS LAI values are shown in Fig. 6. The mean differences between the MODIS and the default LAI values are presented in Fig. 7. The largest positive differences occur in areas with “arable land”, where the mean default LAI values are lower than the MODIS LAI values. The largest negative deviations occur in areas with forest, especially for “coniferous forest”. The seasonal variations of the MODIS and the default LAI values are shown in Fig. 8. The default LAI of “grass” and “deciduous forest” seems to fit the yearly variation of the MODIS LAI quite well. We matched the MODIS LAI with the locations of the FLUXNET sites to take a closer look at the pattern for different land use classes. Figures 9 and S3 show the seasonal variation of the MODIS LAI at FLUXNET sites with different land use classifications. The LAI of the grassland sites seems to vary the most, which corresponds to the large inter-site differences in vegetation cover as explained in Sect. 4.3. For the cropland sites, we can recognize the growing season and the apparent harvest, where the LAI values drop again. Of all the different land use classes, deciduous broadleaf forest sites reach the highest LAI values in the growing season. There is less variation in the LAI for evergreen needleleaf forest sites. However, the wintertime LAI values seem to be unrealistically low.

4.5 Implications for modelled N_r deposition fields

In the following section, the impact of the updated LAI and z₀ 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 of the updated LAI and z₀ 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₀ values. The second run, named “LE_LAI”, uses updated LAI values and the default z₀ values. The third run, named “LE_z0”, uses updated z₀ values and the default LAI values. The fourth run, named “LE_z0+LAI”, uses both updated LAI and z₀ values. From now on, we will refer to the outputs of these different runs using the abovementioned abbreviations.

4.5.1 Effect on total N_r deposition

Figure S4 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 corresponds to the relative increase in LAI values here. The inclusion of the updated z₀ values led 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 %).

4.5.2 Effect on wet and dry N_r deposition

We examined the direct effect of the updated LAI and z₀ values on the modelled dry N_r deposition, as well as the related indirect effect in modelled wet N_r deposition. Figure 10 shows the dry and wet N_r deposition in kg N ha⁻¹ in 2014, modelled with the updated LAI and z₀ 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 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 %.

4.5.3 Effect on reduced and oxidized N_r deposition

We split up the total N_r deposition into NH_x (NH₃ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$) and NO_y (NO and NO₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and HNO₃) deposition to look at the effect of the updated LAI and z₀ input maps on the deposition of reduced and oxidized forms of N_r, respectively. Figure 12 shows the modelled NH_x and NO_y deposition in kg N ha⁻¹ in 2014, including the updated LAI and z₀ 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₀ values have a larger impact on the NH_x deposition than on the NO_y deposition. The implementation of the updated z₀ values has led to a decrease in NH_x 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.

Figure 11 The relative change in total dry (a, b, c) and wet (d, e, f) N_r deposition in 2014 for the different model runs relative to the default LOTOS-EUROS run. The first row indicates the changes related to the implementation of the updated z₀ values. The second row indicates the changes related to the implementation of the MODIS LAI values. The third row shows the combined effect of both these updates. (Please note the different scales).

Figure 12 NH_x deposition (a) and NO_y deposition (b) in kg N ha⁻¹ in 2014.

Table 7 Relative change (%) in total N_r deposition with respect to the default run over Germany, the Netherlands and Belgium in 2014 per land use class.

Figure 13 The relative change in total NH_x (a, d) and NO_y (c, f) deposition in 2014 for the different model runs relative to the default LOTOS-EUROS run. The first row indicates the changes related to the implementation of the updated z₀ values. The second row indicates the changes related to the implementation of the MODIS LAI values. The third row shows the combined effect of both these updates.

4.5.4 Effect on land-use-specific fluxes

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

4.6 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₀ values on NH₃ and NO₂ surface concentrations. Figure 14 shows the updated modelled NH₃ and NO₂ surface concentrations in 2014. The dots on top of the figures represent the stations that are used for validation and their observed mean NH₃ and NO₂ surface concentrations. Figure 15 shows the relative change in yearly mean NH₃ and NO₂ surface concentrations in 2014 of the different runs with respect to the default run of LOTOS-EUROS.

Figure 14 The yearly mean NH₃ (a) and NO₂ (b) surface concentrations in µg m⁻³ in 2014 and the corresponding mean surface concentrations measured at the in situ stations.

Figure 15 The relative change (%) in mean NH₃ (a, b, c) and NO₂ (d, e, f) surface concentration in 2014 for the different model runs relative to the default LOTOS-EUROS run. The first row indicates the changes related to the implementation of the updated z₀ values. The second row indicates the changes related to the implementation of the MODIS LAI values. The third row shows the combined effect of both these updates.

Figure 16 The mean ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (a) and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (b) wet deposition in kg N ha⁻¹ month⁻¹ in 2014. The mean observed wet deposition observed at the stations is plotted on top.

The first column represents the changes in NH₃ and NO₂ surface concentrations due to the updated z₀ values. The NH₃ surface concentration in the Netherlands has increased by up to ∼8 %. The NH₃ surface concentration in almost all of Germany has decreased by up to ∼10 %. The changes in the NO₂ surface concentration are less distinct and changed approximately minus to plus 4 % in most areas. The middle column represents the changes in NH₃ and NO₂ surface concentrations due to the inclusion of the MODIS LAI only. The NH₃ surface concentration has increased with up to ∼10 % in the Netherlands, western Belgium and north-western and southern Germany. The NH₃ surface concentration has decreased in eastern Belgium and central Germany. The NO₂ surface concentrations have decreased with up to ∼6 % in almost all of the modelled area.

4.7 Comparison to in situ measurements

To analyse the effect of the updated LAI and z₀ 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 ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ wet deposition and NH₃ and NO₂ surface concentrations measurements instead. The distribution of the wet deposition stations is shown in Fig. 16, as well as the modelled mean ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (left) and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (right) wet deposition in 2014. The locations of the stations that measure the NH₃ and NO₂ surface concentrations are shown in Fig. 14.

The relationships between the modelled and observed fields are evaluated using the Pearson 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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ wet deposition, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition and the monthly mean NH₃ and NO₂ surface concentrations in 2014. Table S1 shows the same statistics but computed per DEPAC land use class. Overall, we do not observe large changes in the shown measures due to the inclusion of the updated LAI and z₀ values on a yearly basis. The model underestimates ${\mathrm{NO}}_{\mathrm{3}}^{-}$ wet deposition, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to a lesser extent, too. The modelled NH₃ surface concentrations are in general higher than observed concentrations. The opposite is true for NO₂; 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₀ and LAI values.

Figure S5 shows the monthly mean ${\mathrm{NO}}_{\mathrm{3}}^{-}$ wet deposition, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition, NH₃ surface concentration and NO₂ surface concentrations for the different model runs and the mean of the corresponding in situ observations. The relative changes with respect to the default model run are shown in the bottom figures. For ${\mathrm{NH}}_{\mathrm{4}}^{+}$, the model captures the observed pattern quite well, although the mean spring peak has slightly shifted. The model captures the monthly variation of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ well in general, too. There appears to be an underestimation during the winter, especially in December. The observed NH₃ surface concentrations are lower than the modelled concentrations at the beginning of spring and higher during summer. The measured NO₂ surface concentrations are continuously higher than the modelled values. A potential reason for this might be the spread of the NO₂ stations. Unlike NH₃, NO₂ is not only measured in nature areas but all types of locations. Even the selected background stations may therefore be located relatively closer to emission sources, leading to higher observed NO₂ surface concentrations. The changes due to the inclusion of either the MODIS LAI or the updated z₀ values in our model are limited.

Both Tables S1 and S2, and Fig. 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₀ 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 improvements with the use of the currently available in situ measurements.

5 Discussion

This paper aimed to improve the surface characterization of LOTOS-EUROS through the inclusion of satellite-derived LAI and roughness length (z₀) values. We used empirical functions to derive z₀ values for different land use classes. The updated z₀ values are compared to literature values, showing a good agreement in general. We also compared our z₀ values to z₀ values computed from FLUXNET sites. The z₀ values for forests seemed to match well, but the z₀ values for short vegetation seem to be overestimated for crops and underestimated for grassland and wetland sites. The differences for short vegetation types can be partially explained by the large inter-site variability in vegetation types within each classified land use (e.g. reeds versus short grass). The equation for short vegetation used here seems to work best for short grasslands. For our current study area, this does not pose a problem, since most grasslands in Germany, Belgium and the Netherlands are managed and grazed upon. We found an improved RMSD value of 0.60, compared to RMSD of 0.76 with default z₀ values. Even though there is an offset between the satellite-derived and computed FLUXNET z₀ values for crops, the seasonal pattern seemed to match well. The offset can be explained by the absence of low NDVI (<0.4) values.

The z₀ is closely related to the geometric features and distributions of the roughness elements in a certain area. In reality, it is not only dependent on vegetation properties but, for instance, also the topography of the area. The updated z₀ 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₀ 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₀ values on the deposition fields is limited. To illustrate this, the relative change in total dry NH₃ deposition due to the updated z₀ values were computed and shown in Fig. 17. We used increasing threshold percentages to sort the NH₃ deposition on a model pixel level per land use type and fraction. Figure 17 shows that the differences in total NH₃ 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₃ deposition. The change thus appears to be less distinct in pixels that have a 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₀ on a model pixel level is thus likely to cause a levelling effect on the current model resolution. The impact of the updated z₀ 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.

Figure 17 The relative difference (%) in total dry NH₃ deposition in 2014 between the default run (LE_default) and the run with the updated z₀ values (LE_z0), sorted by increasing land use fraction.

Figure 18 Critical load exceedances on forests in kg N ha⁻¹ in 2014. The upper figures show the critical on deciduous (a) and coniferous (b) forest, as modelled with the updated z₀ and LAI values. The lower figures show the absolute differences in critical load exceedance on deciduous (c) and coniferous (d) between the new and the default LOTOS-EUROS run.

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 root mean square error (RMSE) was 4.4 m and R² was 0.7 after removal of seven 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., 2020). For all but one site (DE-Hai), the forest canopy heights from GLAS were lower than 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₀ maps.

The MODIS LAI at the FLUXNET sites showed realistic seasonal variations for most land use classifications, except for relatively low wintertime values for evergreen needleleaf 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., 2010; Wang et al., 2004). The MODIS LAI products have been gradually improving with each update; however, these issues still exist in the newer versions of the product. For the most recent version of the MODIS LAI, version 6, Yan et al. (2016) found an overall RMSE of 0.66 and a R² of 0.77 in comparison with LAI ground observation. More recently, using a different approach, Xu et al. (2018) found a slightly higher RMSE of 0.93 and a R² of 0.77. Some studies (e.g. Tian et al., 2004) have reported an underestimation of the MODIS LAI in presence of snow cover, particularly affecting evergreen forests. With only limited amounts of snowfall, the regions in our study did not suffer from this problem. 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 and for these regions 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 distributions of these LAI values are more realistic than those of the default LAI values used in LOTOS-EUROS. The same holds for the updated, time-variant z₀ values. The representation of the growing season is now more realistic due to their dependence on NDVI and LAI values. Moreover, the z₀ values for forests now also have a clear spatial variation, such as a latitudinal gradient with increasing z₀ 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₀ 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₀ values led to a variation of $\sim \pm \mathrm{8}$ %, and the updated LAI values led to variations of $\sim \pm \mathrm{30}$ % in both fields. The dry deposition fields were most sensitive to changes in z₀ and LAI, as these varied from approximately −20 % to +20 % with the updated z₀ 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 subgrid level. To illustrate the potential consequences on a local scale, we computed the critical load exceedances for deciduous and coniferous forest (Fig. 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 up to +2 kg for deciduous forest and even over −3 kg for coniferous forest.

The uncertainties of the LAI and z₀ 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 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₃ emissions is estimated to be at least ±30 %. This is mainly due to the diverse nature of agricultural emission sources, leading to large spatiotemporal variations. The annual NO_x emissions total hold a lower uncertainty, of approximately ±20 % (Kuenen et al., 2014). Emissions at specific locations, especially for NH₃, 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 satellite-derived NH₃ emissions of large point sources are a factor of 2.5 higher than those given in emission N inventories. The surface exchange parameterization is another source of uncertainty. The complexity of the NH₃ surface exchange schemes in CTMs is usually low compared to the current level of process understanding (Flechard et al., 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 of 2–3.

This work has shown that changes in two of the main deposition parameters (LAI, z₀) 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 the 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 of satellite products to derive LAI and z₀ 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₀ estimates presented here can be linked to any land use map and is as such transferable to many different models and geographical areas.