Uncertainty in the determined nitrogen
deposition in the Netherlands
Status report 2023

RIVM report 2022-0085
R. Hoogerbrugge et al.

Uncertainty in the determined nitrogen
deposition in the Netherlands
Status report 2023

RIVM report 2022-0085

RIVM report 2022-0085

Colophon

© RIVM 2023
Parts of this publication may be reproduced, provided acknowledgement
is given to the: National Institute for Public Health and the Environment,
and the title and year of publication are cited.
RIVM attaches a great deal of importance to the accessibility of its
products. However, it is at present not yet possible to provide this
document in a completely accessible form. If a part is not accessible, it
is mentioned as such. Also see www.rivm.nl/en/accessibility
DOI 10.21945/RIVM-2022-0085
R. Hoogerbrugge (author), RIVM
M. Braam (author), RIVM
K. Siteur (author), RIVM
C. Jacobs (author), RIVM
S. Hazelhorst (author), RIVM
G. Stefess (author), RIVM
E. van der Swaluw (author), RIVM
R. Wichink Kruit (author), RIVM
J. Wesseling (author), RIVM
A. van Pul (author), RIVM
Contact:
Ronald Hoogerbrugge
Centre MIL
ronald.hoogerbrugge@rivm.nl
This investigation was performed by order, and for the account, of the
Dutch Ministry of Agriculture, Nature and Food Quality, within the
framework of Project 36: Program Nitrogen & Nature.

Published by:
National Institute for Public Health
and the Environment, RIVM
P.O. Box 1 | 3720 BA Bilthoven
The Netherlands
www.rivm.nl/en

Page 2 of 93

RIVM report 2022-0085

Synopsis

Uncertainty estimates of the nitrogen deposition in the
Netherlands
Status report 2023
Every year, RIVM produces maps that show how much nitrogen is
deposited in the Netherlands. We refer to this as nitrogen deposition.
RIVM determines the nitrogen deposition by combining model
calculations with measurements. These model calculations are necessary
because measurements cannot be performed everywhere in the
Netherlands. The deposition is influenced by the weather, the roughness
of the terrain and characteristics of the vegetation. These factors give
rise to a degree of uncertainty.
In this technical report, RIVM describes how the uncertainties in terms
of deposition levels are derived. In terms of degree, the uncertainties
are similar to earlier estimates in 2004 and 2010, but they are now
better underpinned.
The uncertainty in a calculated deposition level is expressed as the
probability of the calculated level deviating from the actual level. On the
national scale, the probability that the calculated level deviates by
between less than 20% and up to 30% from the actual level is
considerable (95%). For specific areas in the Netherlands, there is a
considerable probability (95%) that the calculated level deviates by
between less than 60% and up to 70% from the actual level. This latter
uncertainty applies when the deposition is calculated for a very small
area, such as a hectare or a square kilometre. The degree of uncertainty
of to the national nitrogen deposition is lower because uncertainties in
the processes affecting deposition average out on the national scale.
The results also show the largest contributors to the uncertainties in the
overall nitrogen deposition. The most important factor in this regard is
the uncertainty in the dry deposition, which entails nitrogen depositing
on soil or vegetation directly from the air. There can be considerable
variation from area to area and there is insufficient information to
quantify this process.
Keywords: nitrogen deposition, uncertainty, dry deposition, OPS, DN,
measuring and modelling, COTAG

Page 3 of 93

RIVM report 2022-0085

Page 4 of 93

RIVM report 2022-0085

Publiekssamenvatting

Onzekerheid in de bepaalde stikstofdepositie in Nederland
Status report 2023
Elk jaar maakt het RIVM kaarten die aangeven hoeveel stikstof in
Nederland op de bodem neerslaat. Dit noemen we stikstofdepositie. Het
RIVM bepaalt de stikstofdepositie door modelberekeningen met
metingen te combineren. Modelberekeningen zijn nodig omdat niet
overal in Nederland kan worden gemeten. Het weer, de ruwheid van het
land en eigenschappen van de begroeiing beïnvloeden de depositie en
zorgen voor een bepaalde onzekerheid.
In dit technische rapport beschrijft het RIVM hoe de onzekerheden in de
depositiewaarden worden bepaald. De grootte van de onzekerheden zijn
vergelijkbaar met de eerdere schattingen uit 2004 en 2010, maar zijn
nu beter onderbouwd.
De onzekerheid in de berekende depositiewaarde drukken we uit in de
kans dat de berekende waarde een bepaalde afwijking heeft van de
werkelijke waarde. Landelijk gezien is de kans groot (95 procent
betrouwbaarheidsinterval) dat de berekende waarde minder dan 30
procent afwijkt van de werkelijke waarde. Op een specifieke locatie in
Nederland is de kans groot (95 procent betrouwbaarheidsinterval) dat
de berekende waarde minder dan 70 procent afwijkt van de werkelijke
waarde. Deze onzekerheid geldt wanneer de depositie wordt berekend
voor een heel klein gebied, zoals een hectare of een vierkante kilometer.
De kleinere onzekerheid in de landelijke stikstofdepositie komt omdat
onzekerheden in processen die invloed hebben op de depositie, op
landelijke schaal uitmiddelen.
De onzekerheden hebben geen betrekking op berekeningen voor
specifieke projecten maar uitsluitend voor de totale depositie. De
resultaten laten ook zien waar de grootste onzekerheden in de totale
stikstofdepositie door komen. De onzekerheid in de droge depositie is
daar het belangrijkste onderdeel van. Bij droge depositie komt stikstof
direct via de lucht op de bodem of in de vegetatie terecht. De
hoeveelheid kan van plek tot plek sterk verschillen en er is weinig
informatie beschikbaar om dit proces te kwantificeren.
Kernwoorden: stikstofdepositie, onzekerheid, droge depositie, OPS, DN,
meten en modelleren, COTAG

Page 5 of 93

RIVM report 2022-0085

Page 6 of 93

RIVM report 2022-0085

Contents

1
1.1
1.1.1
1.1.2
1.1.3
1.1.4
1.2
1.2.1
1.2.2
1.2.3
1.2.4
1.2.5
1.3
1.4
1.4.1
1.4.2
1.4.3
1.5
1.6

Nederlandse samenvatting Onzekerheid in de
stikstofdepositie — 9
Inleiding — 9
Systeem van emissies en deposities — 9
Beschikbare metingen — 10
Modelberekeningen en onzekerheid — 11
Nederlandse samenvatting en rapport — 13
Onzekerheden in de zes componenten van de stikstofdepositie — 13
Methode — 13
Onzekerheden in de natte depositie — 14
Onzekerheden in de droge depositie van ammoniak — 14
Onzekerheden in de droge depositie NO2, NO3- en NH4+ — 16
Conclusies — 17
Totale onzekerheid — 17
Effecten van interacties tussen depositie en concentratie — 18
De negatieve terugkoppeling — 18
Onderzoeksmethode — 18
Resultaten — 19
Vergelijking GCN/GDN met DN — 21
Discussie, conclusie en aanbevelingen — 22

2
2.1
2.2
2.3

Introduction — 25
Aim of the study — 25
Detailed background — 27
Reading Guide — 29

3
3.1
3.2
3.2.1
3.3
3.4
3.5
3.5.1
3.5.2

Uncertainty in nitrogen deposition estimates per compound — 31
Measurements — 31
Model calculations — 33
Post-processing of OPS outputs — 33
Removing systematic model bias — 34
Uncertainty in the model — 35
Uncertainties in dry deposition — 35
Decomposition of the uncertainty in dry deposition fluxes — 35
Estimation of uncertainty in dry deposition velocities of NH4+, NO3-, and
NO2 — 36
Results — 37
Interpretation of difference between modelled and measured dry
deposition fluxes — 41
Conclusion & Discussion — 42

3.6
3.7
3.8
4
4.1
4.2
4.2.1
4.2.2
4.3
4.4
4.5

Uncertainty in total nitrogen deposition — 43
Uncertainty in total nitrogen deposition — 43
Total uncertainty in total nitrogen deposition — 45
Sum of model biases — 45
Estimation of total uncertainty in total nitrogen deposition — 46
Average uncertainty in nitrogen deposition on Natura 2000
habitats — 46
Uncertainty of the national average (total) nitrogen deposition — 47
Conclusions — 50
Page 7 of 93

RIVM report 2022-0085

5
5.1
5.1.1
5.1.2
5.2
5.2.1
5.2.2
5.2.3
5.2.4
5.3
5.3.1
5.3.2
5.3.3
5.3.4
5.3.5
5.4
6

Impact of atmospheric feedback on sensitivity of nitrogen
deposition to deposition velocity — 51
Introduction — 51
Background — 51
Goal — 52
Method — 52
Model simulations: general — 52
Variation of Vd — 52
Land use-dependent variation of Vd at receptor points — 53
Analysis of results — 54
Results and discussion — 54
Deposition velocity Vd — 54
Atmospheric concentration — 56
Deposition — 57
Carry-over effects — 60
Aggregation to larger spatial scales — 60
Conclusions — 63

6.1
6.2
6.3
6.3.1
6.3.2
6.4

Spatial variability of nitrogen deposition: comparing of GCN/GDN
and DN — 65
Introduction — 65
Methods — 65
Results — 66
Comparing of modelled NH3 concentrations with observations — 66
Comparing of DN hexagonal grid and GDN square grid. — 68
Conclusions — 72

7

General discussion and conclusion — 73

8

References — 75
Appendix A Proportional random uncertainty and alternative
models — 79
Appendix B Independence between deposition velocities and
concentrations — 81
Appendix C Uncertainties in dry deposition fluxes and
correlations between fluxes — 82
C.1 Variance of the product of two independent random
variables — 82
C.2 Exact Covariance of Products of Random Variables — 82
C.3 Correlation between residuals in concentrations, dry deposition
velocities and wet deposition fluxes — 83
Appendix D The average of multiple draws of a random
variable — 86
Appendix E — 87
E.1 Assessing Vd of NH3 in an international context — 87
E.2 Interpretation and processing of data in Schrader & Brümmer
(2014) and Vendel et al. (2023) — 88

Page 8 of 93

RIVM report 2022-0085

1

Nederlandse samenvatting
Onzekerheid in de stikstofdepositie

1.1

Inleiding

1.1.1

Systeem van emissies en deposities
Het RIVM brengt elk jaar kaarten uit van de stikstofdepositie in
Nederland. De doelen van deze kaarten zijn beleidsmonitoring en het
monitoren van de stikstofdepositie. Voorbeelden zijn de Grootschalige
Depositiekaarten Nederland, GDN (Hoogerbrugge et al., 2022) en de
Depositiekaarten stikstof Natura 2000-gebieden, DN. De resultaten van
de DN-kaarten worden gebruikt in de stikstofmonitoring voor de Wet
stikstofreductie en natuurverbetering, Wsn (Marra et al., 2022).
De kaarten worden gemaakt op basis van een combinatie van metingen
en modelberekeningen van stikstof.
De modelberekeningen beschrijven de keten van emissie via transport
naar depositie (Figuur 1). Door emissie komen er stikstofverbindingen in
de lucht. Verkeer en industrie stoten voornamelijk stikstofoxiden
(geoxideerd stikstof) uit, en de landbouw voornamelijk ammoniak
(gereduceerd stikstof). De hoeveelheid van een stof per volume-eenheid
heet de concentratie. In de lucht worden de stoffen verder met de wind
mee getransporteerd en vinden er chemische omzettingen plaats.

Figuur 1 Overzicht van de stikstofroutes door de lucht

De stoffen verdwijnen vervolgens uit de lucht als gevolg depositie. Er
zijn twee manieren van depositie (zie de rechterkant in Figuur 1).
Ten eerste kunnen de componenten met het regenwater mee neerslaan
(natte depositie) in de vorm van nitraat (geoxideerd stikstof in
oplossing) en ammonium (ammoniak in oplossing).
Ten tweede kunnen planten, bodem en de waterlaag daarop ze direct uit
de lucht opnemen (droge depositie). Hierbij onderscheiden we
primaire componenten en secundaire componenten. De primaire
componenten zijn direct uitgestoten: ammoniak (NH3) en stikstofoxiden
(NOx = NO + NO2). De secundaire componenten ontstaan door
Page 9 of 93

RIVM report 2022-0085

chemische reacties in de atmosfeer (bijvoorbeeld ammoniumsulfaat,
ammoniumnitraat).
Dit rapport gaat dan ook over zes componenten waaruit de totale
stikstofdepositie bestaat:
𝐷𝐷N tot = 𝐷𝐷NO3 nat + 𝐷𝐷NOx primair droog + 𝐷𝐷NOx secundair droog +
𝐷𝐷NH4 nat + 𝐷𝐷NH3 primair droog + 𝐷𝐷NH4 secundair droog

Vgl. 1

Figuur 2 Gemiddelde berekende Nederlandse depositiebijdrage per component
over de periode 2005-2021 voor geoxideerd (links) en gereduceerd stikstof
(rechts). Voor 𝐷𝐷𝑁𝑁𝑁𝑁𝑥𝑥 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , 𝐷𝐷𝑁𝑁𝑁𝑁3 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , 𝐷𝐷𝑁𝑁𝑁𝑁4 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 (gestreept) zijn
alleen concentratiemetingen beschikbaar. Voor 𝐷𝐷𝑁𝑁𝑁𝑁3 𝑛𝑛𝑛𝑛𝑛𝑛 , 𝐷𝐷𝑁𝑁𝑁𝑁3 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , 𝐷𝐷𝑁𝑁𝑁𝑁4 𝑛𝑛𝑛𝑛𝑛𝑛
(gevuld en zwart omkaderd) zijn concentratie- én depositiemetingen
beschikbaar.

In Figuur 2 staat hoeveel iedere component bijdraagt aan de totale
Nederlandse depositie. Verder is aangegeven of van een component
alléén concentratiemetingen beschikbaar zijn, of zowel concentratie- als
depositiemetingen.
1.1.2

Beschikbare metingen
Een overzicht van alle relevante stikstofmetingen staat in het rapport
over de meetstrategie stikstof (Wichink Kruit et al., 2021). Van alle
stikstofcomponenten zijn jaarlijkse concentratiemetingen beschikbaar.
Het meten van natte depositie is relatief eenvoudig. Natte depositie
meten we door eerst het regenwater op te vangen, en dan te meten
hoeveel er van de stikstofcomponent in het regenwater zit.
Het meten van droge depositie van gasvormige componenten en
fijnstofdeeltjes is moeilijker. Figuur 2 laat zien dat de droge depositie
van ammoniak de grootste bijdrage heeft in Nederland.
Het RIVM monitort op zes plekken droge depositie van ammoniak met
een meetsysteem dat we het COTAG-systeem noemen (Rutledge-Jonker
et al., 2023). In Bargerveen (2012), de Oostelijke Vechtplassen (2014),
de Hoge Veluwe (2017), het Noordhollands Duinreservaat (2021),
Veenkampen (2022) en Haaksbergerveen (2023). Er waren tijdens de
analyse van dit onderzoek alleen van de eerste drie COTAG-systemen
jaarcijfers beschikbaar (Rutledge-Jonker et al., 2023). In de analyses
van dit rapport gebruiken we alleen die metingen. Van ammoniak
hebben we dus concentratie- én enkele depositiemetingen.
Page 10 of 93

RIVM report 2022-0085

De droge depositie van stikstofoxides en fijnstofdeeltjes wordt nu (nog)
niet gemeten. Voor deze componenten zijn dus alleen de
concentratiemetingen beschikbaar.
1.1.3
1.1.3.1

Modelberekeningen en onzekerheid
Onzekerheden
Het is niet mogelijk om overal in Nederland te meten. Om toch de
situatie voor heel Nederland in kaart te brengen, gebruikt het RIVM
rekenmodellen. Deze modellen worden ook gebruikt om te berekenen
hoe de concentratie en depositie van stikstof zich in de toekomst kunnen
ontwikkelen. Deze modelberekeningen hebben een bepaalde mate van
onzekerheid.
Het bestuderen van de onzekerheid in de stikstofdepositie is lastiger dan
van de concentraties in de lucht. Dit komt doordat depositie moeilijker
te meten is. Maar ook het modelleren van het depositieproces van
stikstof is complex. Het model moet alle zes de componenten en de
bijhorende depositiepaden meenemen en berekenen.
De onzekerheid in de gemodelleerde totale stikstofdepositie op
1 x 1 km2 is eerder geschat op 70% (Velders et al., 2010). Deze waarde
was gebaseerd op een analyse uit 2004 (van Jaarsveld, 2004). In dit
rapport geven we een update van deze waarde, gebaseerd op een
uitbreiding van de eerder gevolgde methode. Schattingen van de
onzekerheden in concentraties en de natte depositie zijn beschikbaar,
bijvoorbeeld in de EMEP-statusrapporten (EMEP, 2022). Dit geldt niet
voor schattingen van de onzekerheid in de totale stikstofdepositie.
Walker et al. (2019) schrijven dat er geen generieke schattingen
beschikbaar zijn van de onzekerheid in gemodelleerde stikstofdeposities.

1.1.3.2

Modelprocessen
Zoals we eerder schreven, hangt de depositie af van de concentratie.
Voor de natte depositie neemt het model in de berekening ook de
neerslagduur en -intensiteit mee.
De hoeveelheid droge depositie is een vermenigvuldiging van de
luchtconcentratie van de stikstofcomponent met de depositiesnelheid
van die component. De depositiesnelheid is een maat voor het gemak
waarmee stoffen door de lucht naar het aardoppervlak bewegen en
daarna door bodem en planten worden opgenomen. Dit is voor elke
stikstofcomponent anders. Belangrijke factoren die depositiesnelheid
bepalen, zijn het weer, landgebruik, ruwheid, de bodemsoort en
vegetatiesoort.
Bij de droge depositie van ammoniak op vegetatie is ook de hoeveelheid
ammoniak die in de vegetatie zit van belang. Modellen houden hier
rekening mee door het zogenoemde compensatiepunt te bepalen. Deze
waarde hangt af van de plantensoort en bijvoorbeeld de temperatuur.
Als het compensatiepunt hoger is dan de concentratie in de lucht is er
geen depositie maar emissie.
We gaven al aan hoe de hoeveelheid concentratie een effect heeft op de
depositie, maar het omgekeerde geldt ook. De hoeveelheid depositie
heeft invloed op de concentratie. Depositie draagt bij aan de
verwijdering van een component uit de atmosfeer, waardoor de
Page 11 of 93

RIVM report 2022-0085

concentratie daalt. Deze interactie tussen concentratie en depositie
wordt ook wel een negatieve terugkoppeling genoemd. Dit bestuderen
we in meer detail in paragraaf 1.4 en Hoofdstuk 5.
1.1.3.3

Het OPS-model, GCN/GDN en DN
Het RIVM berekent de stikstofdepositie met het Operationele Prioritaire
Stoffen model (OPS-model; Sauter et al., 2023; van Jaarsveld, 2004).
Dit model is een rekenprogramma om de verspreiding van
verontreinigende stoffen in de lucht en de depositie te modelleren. OPS
is de basis voor de GCN/GDN-kaarten (Hoogerbrugge et al., 2022) en
voor de DN-kaarten (https://monitor.aerius.nl/).
Voor beide modelsystemen leggen we een rooster (‘grid’) van vakjes
over Nederland en berekenen we voor elk vakje de concentratie en
depositie van de verschillende componenten. De vorm en schaal
(resolutie) van de vakjes verschilt tussen de modellen. GCN/GDN rekent
op vierkanten van 1 × 1 km2. DN rekent op hexagonen (zeshoeken) met
een oppervlak van 1 ha. GCN/GDN rekent voor heel Nederland, en DN
rekent specifiek voor de stikstofgevoelige natuur in Natura 2000gebieden.

1.1.3.4

De modelkalibratie
In GCN/GDN en DN kalibreren we de modeluitkomsten met metingen.
Dit wordt gedaan om een systematische afwijking tussen het model en
de metingen te corrigeren. De aanpak van de kalibratie verschilt per
component:
• De natte depositie kalibreren we aan de hand van de natte
depositiemetingen.
• Voor de droge depositie van ammoniak kalibreren we aan de
hand van alleen de concentratiemetingen (Hoogerbrugge et al.,
2022; Velders et al., 2010). Dit komt doordat er onvoldoende
droge depositiemetingen beschikbaar zijn en omdat de metingen
een lokale weergave zijn en dus minder geschikt voor landelijke
interpolatie..
Met andere woorden: de kalibratie van de droge depositie is op
één deel van het proces uitgevoerd, namelijk alleen de
concentratie en niet de depositiesnelheid. Dit is niet optimaal,
want een belangrijke processtap kan dus niet worden
gekalibreerd.
• De droge depositie van NO2 en secundaire componenten worden
niet gekalibreerd in GCN/GDN en DN.
In deze Nederlandse samenvatting rapporteren we de gekalibreerde
modelberekeningen. Daarom zijn de onzekerheden ook representatief
voor de nauwkeurigheid van de gepubliceerde deposities in GCN/GDN en
DN. De berekende onzekerheden betreffen de onzekerheden in de totale
depositie en hebben geen betrekking op de onzekerheid van een
projectberekening voor een individuele bron.
In het rapport tonen we ook de ongekalibreerde resultaten, omdat die
informatie relevant is voor het mogelijke verbeteringen in het OPSmodel.

Page 12 of 93

RIVM report 2022-0085

1.1.4

Nederlandse samenvatting en rapport
We bespreken in dit rapport de volgende zaken:
1. De onzekerheden in de zes componenten en de totale depositie.
Deze onzekerheden bepalen we door de modelwaarden te
vergelijken met metingen (Paragraaf 1.2 en Hoofdstuk 3).
2. De totale onzekerheid van de zes componenten voor Nederland
en voor de Nederlandse Natura 2000-gebieden (Paragraaf 1.3 en
Hoofdstuk 4).
3. De negatieve terugkoppeling tussen de droge depositie en de
concentratie (Paragraaf 1.4 en Hoofdstuk 5)
4. Een vergelijking tussen GCN/GDN en DN (Paragraaf 1.5 en
Hoofdstuk 6)
5. Ten slotte de discussie, conclusie en aanbevelingen (Paragraaf
1.6 en Hoofdstuk 7). Dit rapport is een ‘statusrapport’. Het geeft
de huidige stand van zaken weer over de onzekerheden in de
modellering van de stikstofdepositie. In dit laatste deel
bespreken we de beperkingen en verbeterpunten.

1.2

Onzekerheden in de zes componenten van de stikstofdepositie

1.2.1

Methode
De eerste stap is het bepalen van de onzekerheid voor de zes
stikstofcomponenten. We doen dit door jaarcijfers van de OPSberekeningen te vergelijken met de beschikbare metingen. We nemen
data van alle beschikbare jaren mee in de periode van 2005 tot 2021.
Zoals we eerder aangaven, rapporteren we in de Nederlandse
samenvatting alleen de gekalibreerde modelberekeningen. Voor de
meeste componenten is de modelkalibratie ook gebaseerd op metingen
(zie vorige paragraaf, Hoogerbrugge et al., 2022). Met andere woorden,
we gebruiken de meting dus zowel voor de modelkalibratie als voor de
modelvalidatie. Om ervoor te zorgen dat we de validatie onafhankelijk
van de kalibratie uitvoeren, gebruiken we de ‘leave one out’-methode.
Dit betekent dat we tijdens de validatie voor elk meetpunt een
kalibratiekaart maken waar dat specifieke meetpunt niet in zit. Na
kalibratie is het gemiddelde verschil tussen modelberekening en meting
weggenomen en liggen de punten rond de 1:1-lijn. De onzekerheid die
wij hier beschouwen, is dus de afwijking ten opzichte van deze lijn.
De onzekerheden van de zes componenten hieronder rapporteren we in
2𝜎𝜎. Dit komt voor een normale verdeling ongeveer overeen met het
95% betrouwbaarheidsinterval. Stel, de 2𝜎𝜎-onzekerheidswaarde is 34%
en de meting is 400 mol/ha/y (zie de resultaten voor natte depositie in
de volgende paragraaf). Dan is er 2,5% kans dat de modelwaarde
eigenlijk 34% lager is, dus 256 mol/ha/y. En er is ook 2,5% kans dat de
modelwaarde 34% hoger is, dus 536 mol/ha/y. We bepalen de
onzekerheden door de berekeningen te vergelijken met de metingen. In
de modelonzekerheid is dus ook de meetonzekerheid meegenomen. De
impact van de meetonzekerheid is in dit rapport ook geanalyseerd.

Page 13 of 93

RIVM report 2022-0085

Figuur 3 Vergelijking tussen gemeten en berekende natte depositie van NH4+ (a)
en NO3- (b).

1.2.2

Onzekerheden in de natte depositie
Voor de natte depositie zijn metingen beschikbaar. Dus hier kunnen we
metingen en modelberekeningen direct met elkaar vergelijken. In
Figuur 3 staat deze relatie voor ammonium (NH4+) en geoxideerd
stikstof (nitraat, NO3-).
Uit de spreiding rondom de 1:1-lijn volgt dat de relatieve
modelonzekerheid (2𝜎𝜎) 34% is voor ammonium, en 26% voor nitraat.

1.2.3

Onzekerheden in de droge depositie van ammoniak
Van de zes componenten is voor Nederland de gasvormige droge
depositie van ammoniak de grootste. Dit is ook de enige component
waarvan momenteel droge depositiemetingen beschikbaar zijn.
Daarnaast zijn er voor ammoniak concentratiemetingen beschikbaar.
Hierdoor kunnen we ook een schatting maken van de onzekerheid in de
depositiesnelheid van ammoniak.

Page 14 of 93

RIVM report 2022-0085

Figuur 4a Vergelijking tussen de gemeten en berekende droge depositie van
NH3. De kleuren geven de verschillende locaties aan: oranje is BAR, Bargerveen,
blauw is OVP, Oostelijke Vechtplassen en bruin VEL, de Hoge Veluwe.
Figuur 4b Vergelijking tussen de gemeten en berekende NH3-concentraties. De
kleuren geven de aantal metingen in de hexagoon weer: hoe donkerder hoe
meer metingen.

De vergelijkingen tussen de modelberekeningen en de metingen voor de
droge depositie en concentratie van ammoniak staan in Figuur 4.
De vergelijking voor droge depositie van ammoniak tussen de COTAGmetingen en het OPS-model laat een grote spreiding zien (Figuur 4a).
Hierdoor is de bijhorende onzekerheid ook groot: ongeveer 120% (2𝜎𝜎).
In de meeste gevallen zijn de meetwaarden hoger dan de modelwaarde.
Ze liggen onder de 1:1-lijn. Het is niet mogelijk om een modelkalibratie
voor de droge depositie uit te voeren, aangezien het hier om maar drie
locaties gaat.
Uit de concentratie- en depositiemetingen hebben we ook
depositiesnelheid van ammoniak afgeleid. Deze waarden hebben we ook
vergeleken met eerdere afgeleide depositiesnelheden (Appendix E). Dit
zijn vaak buitenlandse metingen. Deze literatuurwaarden laten een
grote spreiding zien, waar de recente Nederlandse metingen bij
aansluiten. Deze grote spreiding is een tweede reden dat het niet
verstandig is om op basis van een kleine set metingen generieke
kalibraties toe te passen.
Van de ammoniakconcentratie zijn grote sets aan metingen beschikbaar.
Deze vergelijking staat in Figuur 4b. Uit de vergelijking tussen de
berekende en gemeten NH3-concentraties volgt een onzekerheid van
45% (2𝜎𝜎).
OPS berekent de droge depositie door de depositiesnelheid met de
concentratie te vermenigvuldigen. Door deze twee onzekerheden te
combineren kunnen we indirect een onzekerheid afleiden van de
effectieve depositiesnelheid voor ammoniak. Die is 113% (2𝜎𝜎).

Page 15 of 93

RIVM report 2022-0085

1.2.4

Onzekerheden in de droge depositie NO2, NO3- en NH4+
Voor stikstofdioxide en de fijnstofdeeltjes die stikstof bevatten (NO3- en
NH4+) zijn geen recente Nederlandse droge depositiemetingen
beschikbaar. Hier vergelijken we alleen de berekende en gemeten
concentraties met elkaar (Figuur 5). De onzekerheid in de
modelconcentratie is 21% voor NO2, 45% voor NO3- en 72% voor NH4+
(2𝜎𝜎).

Figuur 5a Vergelijking tussen de gemeten en berekende NH4+-concentraties.
Figuur 5b Vergelijking tussen de gemeten en berekende NO2 -concentraties. De
kleuren geven de aantal metingen in de hexagoon weer: hoe donkerder hoe
meer metingen.c: Vergelijking tussen de gemeten en berekende NO3-concentraties.

Voor de schatting van de onzekerheid in de droge depositie is er ook een
schatting in de onzekerheid van de depositiesnelheid nodig. We nemen
de berekende onzekerheid in de depositiesnelheid van ammoniak, en
trekken daar de onzekerheid door het compensatiepunt van af (zie
§3.5.2). De reden hiervoor is dat NO2 in tegenstelling tot NH3 geen
compensatiepunt heeft. Verder hangt de bepaling van de
compensatiepunt van verschillende parameters af en is deze daardoor
Page 16 of 93

RIVM report 2022-0085

onzeker. Uiteindelijk komt de onzekerheid van de depositiesnelheid van
NO2 uit op 91% (2𝜎𝜎). En de onzekerheid in de droge depositie van
stikstofoxides op 94% (2𝜎𝜎).
Voor de onzekerheid in droge depositiesnelheden van fijnstofdeeltjes
gebruiken we dezelfde waarde van 91% (2𝜎𝜎). Hierdoor wordt de
onzekerheid in droge depositie van NH4+ 121% (2𝜎𝜎) en van NO3- 104%
(2𝜎𝜎).
1.2.5

Conclusies
Deze analyse geeft aan dat de onzekerheid in de droge depositiesnelheid
de grootste bijdrage levert aan de onzekerheid in de totale depositie.

1.3

Totale onzekerheid
Nu de onzekerheden van alle zes componenten bekend zijn, kunnen we
ze combineren tot een totale onzekerheid.
In de onzekerheidsanalyse van de totale depositie kunnen we niet
aannemen dat de onzekerheden in de zes componenten onafhankelijk
van elkaar zijn. Sommige componenten correleren met elkaar. Zo wordt
natte depositie beïnvloed door de regenduur en regenintensiteit, en
daardoor zal de natte depositie van ammonium en nitraat voor een deel
gecorreleerd zijn.
We maken daarom een correlatietabel (Table 5 in Paragraaf 4.1), waar
voor alle componenten onderling een correlatie is weergegeven.
De onderlinge correlatie tussen de natte depositiecomponenten bepalen
we uit de metingen, want dit zijn gecombineerde metingen. Voor de
andere componenten schatten we de correlatie met behulp van
modelberekeningen (zie voor meer detail Paragraaf 4.1).
Deze tabel gebruiken we bij het optellen van de onzekerheden, om zo
tot een totale onzekerheid in de stikstofdepositie te komen. De totale
onzekerheid is dan een combinatie van de onzekerheden per component
en de onzekerheid door de correlatie tussen de componenten.
Dit resulteert in een gemiddelde onzekerheid in de totale
stikstofdepositie van ongeveer 60-70% (2𝜎𝜎) voor een willekeurige
locatie in Nederlands natuurgebied. Dit betekent 900-1100 mol/ha/j.
Deze 60% à 70% komt goed overeen met de eerdere schatting (Velders
et al., 2010). Het voordeel ten opzichte van de eerdere
onzekerheidsschattingen is dat de definitie van de onzekerheidsmaat nu
veel scherper is. Merk op dat dit de onzekerheid is van de waarde zelf,
dus in absolute zin. Er is 5% kans dat de gerapporteerde modelwaarde
60-70% lager of hoger is. Maar in veel gevallen worden modelwaarden
relatief met elkaar vergeleken. Voorbeelden hiervan zijn vergelijking van
verschillende jaren in trendanalyses, of het vergelijken van het wel of
niet meenemen van bepaalde bronnen. In deze situaties vallen veel van
de bovenstaande onzekerheden tegen elkaar weg. Hierdoor zal de
onzekerheid kleiner worden. Dit wordt in een volgende rapportage
verder uitgewerkt.

Page 17 of 93

RIVM report 2022-0085

1.4

Effecten van interacties tussen depositie en concentratie

1.4.1

De negatieve terugkoppeling
In de vorige paragrafen bleek dat de onzekerheid in de droge depositie
wordt gedomineerd door de onzekerheid in de depositiesnelheid. In
werkelijkheid is de relatie tussen de onzekerheid in de depositie en de
depositiesnelheid complexer.
Laten we de vergelijking van de depositie erbij pakken. De droge
depositie (𝐷𝐷droog ) is de concentratie (𝐶𝐶) vermenigvuldigd met de
depositiesnelheid (𝑉𝑉d ):
Vgl. 2

Een lokale verandering in de depositie zal de lokale concentratie in de
atmosfeer nauwelijks beïnvloeden. Maar dit wordt anders op een grotere
ruimtelijke schaal, vanaf enkele kilometers.
Stel, de depositiesnelheid is te hoog gemodelleerd. Volg de
genummerde pijlen in Vgl. 2. Merk op dat het ook omgekeerd geldt.
Dan is er daardoor ook een hogere depositie.
Bij transport over langere afstanden zal deze hogere depositie
leiden tot merkbaar lagere concentraties. Want bodem en
vegetatie nemen meer stoffen uit de lucht op.
Daardoor zal de depositie weer enigszins afnemen. Want er zitten
minder stoffen in de lucht die opgenomen kunnen worden.
De interactie tussen depositie, depositiesnelheid en concentratie zorgt
dus voor een negatieve terugkoppeling. De initiële toename van de
depositie is uiteindelijk minder groot, dankzij de extra afname van de
concentratie tijdens het transport. Dit betekent een demping van de fout
of de onzekerheid in de depositie op grotere schaal.
1.4.2

Onderzoeksmethode
In Hoofdstuk 5 onderzoeken we in hoeverre dit effect de onzekerheden
in de depositie beïnvloedt. Dat doen we door verschillende berekeningen
(‘runs’) met OPS uit te voeren.
In elke run gebruiken we een andere depositiesnelheid voor alle
graslanden in Nederland. We hebben de depositiesnelheid beïnvloed
door de oppervlakteweerstand van de graslanden aan te passen. 1 In de
Nederlandse samenvatting laten we alleen de resultaten zien van de run
waar de oppervlakteweerstand 20% is van normaal en de run waar die
500% is van normaal. In Hoofdstuk 3 staan de resultaten van de
overige runs.
Daarna berekenen we op ruim driehonderd verschillende locaties de
lokale depositiesnelheid en de depositie ter plaatse. We onderzoeken
dan per run wat de invloed is van de toename of afname van de
oppervlakteweerstand van gras op de lokale depositiesnelheid en

1
Zoals in 1.1.3 is beschreven, hangt de depositiesnelheid af van verschillende factoren, zoals het weer, het
gewas en de bodem. OPS berekent de depositiesnelheid met het weerstandsmodel. Eén van de drie parameters
van dit model is de zogenaamde oppervlakteweerstand. Als de oppervlakteweerstand laag is, nemen de bodem
(of het water) en de vegetatie gemakkelijk de bepaalde componenten op, en is de depositiesnelheid hoog. Zie
ook Hoofdstuk 5.

Page 18 of 93

RIVM report 2022-0085

depositie. Ook kijken we naar het effect op de totale depositie in heel
Nederland.
We verwachten dat het aanpassen van de depositiesnelheid via de
oppervlakteweerstand van graslanden een effect heeft op de totale
depositie in heel Nederland. Dit komt doordat bijna 40% van het
Nederlandse oppervlak grasland is. De depositieverandering heeft door
transport met de wind niet alleen boven grasland een effect op de
concentratie in de lucht, maar ook boven ander landgebruik.

Figuur 6 Het relatieve effect van de verandering van de oppervlakteweerstand
op de lokale depositiesnelheid (horizontale as) en op de berekende totale
depositie (verticale as) van NHx (links) en NOy (rechts).
De totale depositie is de som van de droge en natte depositie van alle
componenten, zoals berekend met OPS (𝐷𝐷𝑁𝑁𝑁𝑁𝑥𝑥 𝑡𝑡𝑡𝑡𝑡𝑡 = 𝐷𝐷𝑁𝑁𝑁𝑁4 𝑛𝑛𝑛𝑛𝑛𝑛 + 𝐷𝐷𝑁𝑁𝑁𝑁 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 +
𝐷𝐷𝑁𝑁𝑁𝑁4 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , en 𝐷𝐷𝑁𝑁𝑁𝑁𝑦𝑦 𝑡𝑡𝑡𝑡𝑡𝑡 = 𝐷𝐷

𝑁𝑁𝑁𝑁3 𝑛𝑛𝑛𝑛𝑛𝑛

3

+ 𝐷𝐷𝑁𝑁𝑁𝑁𝑥𝑥 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 + 𝐷𝐷𝑁𝑁𝑁𝑁3 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 ).

Bij oranje run is de oppervlakteweerstand van gras 20% van de normale
waarden, en bij de blauwe run 500%.
Veranderingen worden weergegeven als de verhouding met het resultaat van de
normale run. De + tekens geven een uit deze berekeningen geschat gemiddeld
effect voor Nederland, na weging met landgebruik.

1.4.3

Resultaten
De OPS-resultaten laten diverse vormen van compensatie zien. Figuur 6
geeft weer hoe een verandering in depositiesnelheid de lokale en
landelijk gemiddelde depositie beïnvloedt. De cirkels zijn de ruim
driehonderd locaties waar we de lokale depositiesnelheid en depositie
bepalen. Op de horizontale as staat de verhouding in de lokale
depositiesnelheid van de nieuwe run ten opzichte van normaal. Dit geeft
dus de lokale invloed van de aanpassing van de oppervlakteweerstand
van gras op de depositiesnelheid weer. Deze verhouding zal afhangen
van de hoeveelheid gras direct rondom een specifieke locatie. Op de
verticale as staat de verhouding van de totale depositie ten opzichte van
normaal, ofwel het effect op de depositie.
De figuur laat de resultaten zien van twee runs op alle locaties. Bij de
oranje run is de oppervlakteweerstand van gras 20% van de normale
waarden (een factor 5 lager). Hierdoor is de depositiesnelheid boven
gras hoger dan normaal (in de meeste gevallen is de verhouding
𝑉𝑉𝑑𝑑 > 1.0). De invloed op 𝑉𝑉𝑑𝑑 verschilt per locatie omdat de hoeveelheid

Page 19 of 93

RIVM report 2022-0085

grasland in de omgeving van elke locatie verschilt. Bij de blauwe run is
de oppervlakteweerstand 500% (een factor 5) groter dan normaal. We
zien nu het omgekeerde, de depositiesnelheid is lager dan normaal.
De invloed op de depositie is lager dan de verandering van de
depositiesnelheid: de verticale bandbreedte is kleiner dan horizontale.
Zie bijvoorbeeld de oranje 20%-run van NOy: het effect op 𝑉𝑉𝑑𝑑 NOy ligt
tussen 1.0 en 2.0 en voor 𝐷𝐷tot NOy tussen 0.9 en 1.5. Dit komt door de
invloed van de verandering op de concentratie in de lucht. De hogere
depositiesnelheid van grasland in heel Nederland gaat gepaard met
lagere concentraties in de lucht en omgekeerd. De droge depositie op de
individuele locaties wordt daardoor iets lager dan je zou verwachten op
basis van de verandering in de depositiesnelheid ter plekke.
De concentratieverandering beïnvloedt ook de natte depositie. De lagere
concentratie bij hogere depositiesnelheid leidt tot minder natte
depositie, en omgekeerd. Dit draagt sterk bij aan het dempende effect
op de berekende totale depositie (nat en droog).
In Figuur 6 is ook de invloed op het Nederlands gemiddelde te zien (+).
We zien dat forse veranderingen in de depositiesnelheid nauwelijks
effect hebben op de berekende gemiddelde depositie. Dit komt doordat
er naast de dempende effecten ook concentratieveranderingen boven
gebieden met andere landgebruik (niet-gras) optreden. De bijbehorende
depositieveranderingen zijn in veel gevallen tegengesteld aan die boven
gras. En dus zal dit bij middeling over grote gebieden het effect op de
depositie door de aanpassingen van gras verder compenseren.
De interactie tussen depositie en concentratie op grote schaal zorgt
ervoor dat ook lokaal de onzekerheid in de depositie minder is dan de
onzekerheid in de depositiesnelheid. Voor de stikstofoxides is dit in deze
berekeningen ongeveer een factor twee lager, en voor ammoniak een
factor drie. Na het middelen naar landelijke schaal, heeft een forse
verandering in de depositiesnelheid van gras zelfs nauwelijks effect op
de depositie. De verhouding in de totale gemiddelde depositie blijft in
deze berekeningen vrijwel 1 voor ammoniak. Voor NOx is de relatieve
afwijking in de totale depositie maar ongeveer een kwart van de
gemiddelde relatieve afwijking in de depositiesnelheid.
Dit betekent dat de aangenomen één-op-één foutenvoortplanting in
paragraaf 1.2 en 1.3 in feite schaalafhankelijk is. De onzekerheid neemt
dus af bij het modelleren van grotere gebieden.
Het ligt dan ook voor de hand dat de onzekerheid in de totale
stikstofdepositie op Nederlandse schaal aanzienlijk kleiner is dan de 6070% op lokale schaal.
Volgens de bovenstaande berekeningen zou de fout op Nederlandse
schaal wel eens een factor twee (tot drie) kleiner kunnen zijn. Voorlopig
gebruiken we die factor en komen zo uit op een onzekerheid van 2030% (2𝜎𝜎). Maar een meer definitieve kwantificering van het schaaleffect,
waarbij we ook rekening houden met andere typen afwijkingen en met
meer oppervlaktetypen, vereist nog veel meer verschillende
berekeningen. Dit onderzoek is onderdeel van het Nationaal
Kennisprogramma Stikstof (NKS).

Page 20 of 93

RIVM report 2022-0085

1.5

Vergelijking GCN/GDN met DN
Ten slotte vergelijken we de resultaten van de GCN/GDN-kaarten met
die van de DN-kaarten. Ook hier gebruiken we voor beide modellen de
aan de metingen gekalibreerde modelwaarden. Zoals eerder
beschreven, worden in DN de concentraties en deposities op een
honderdmaal hogere resolutie dan GCN/GDN. Daarom worden in de DN
sommige emissies ook op hogere resolutie gebruikt, zonder de som van
de lokale emissies te veranderen.

Figuur 7 vergelijking tussen berekeningen van GDN op 1 × 1 km2 en DN op
100m2 hexagonen voor de totale stikstofdepositie (som van de zes
componenten). De kleuren geven het aantal waarden binnen een hexagoon in de
grafiek weer: hoe donkerder hoe meer waarden.

Figuur 7 geeft de vergelijking van de depositieberekeningen van DN met
die van GCN/GDN weer. Hierbij hebben we alleen de 1 × 1 km2 gridcellen van GCN/GDN meegenomen, waarbij voor DN minimaal negentig
van de honderd hexagonen beschikbaar zijn. Het gemiddelde van de
DN-waarde is telkens vergeleken met de GDN-waarde. Het
systematische verschil tussen de modellen voor depositie is minder dan
1%. Dit komt doordat hierbij van hetzelfde model, emissie en metingen
uitgegaan is. De verschillen hebben een spreiding van 25 mol/ha
(RMSE) en worden onder andere veroorzaakt door de hogere resolutie
van de DN-emissies. De grootste verschillen treden vooral op als er bij
de meetpunten lokale bronnen zijn.
De variabiliteit van de DN-hexagoondeposities binnen een 1 × 1 km2
grid-cel is aanzienlijk (circa 30% (2𝜎𝜎)).

Page 21 of 93

RIVM report 2022-0085

1.6

Discussie, conclusie en aanbevelingen
Dit statusrapport geeft een overzicht van de beschikbare kennis over
onzekerheden in stikstofdepositieberekeningen van de
monitoringskaarten (GCN/GDN en DN). De onzekerheidsberekening is
gebaseerd op de vergelijking van modelberekeningen met de
beschikbare metingen. De bepaling van de onzekerheid in de depositie is
complexer dan de bepaling van die in concentraties. Dit komt doordat
directe metingen van de droge depositie grotendeels ontbreken.
De stikstofdepositie is opgebouwd uit verschillende componenten. Deze
zijn zowel nat als droog en betreffen zowel ammoniak als stikstofoxides
en een aantal reactieproducten. Gezien de complexiteit is er nog geen
volledige onzekerheidsanalyse bekend. Daarom is dit de eerste aanzet
om tot zo’n schatting te komen. Deze onzekerheid is groter dan de
onzekerheden in bijvoorbeeld de luchtkwaliteitsmodellering.
De belangrijkste conclusies zijn dat:
-de onzekerheid van de totale stikstofdepositie bedraagt ongeveer 6070% (2𝜎𝜎) op een willekeurige locatie in Nederland. Deze schatting is
ongeveer gelijk aan eerdere schattingen. Deze nieuwe schatting is dus
consistent met de informatie die aanwezig was toen er beleidskeuzes
gemaakt werden ten aanzien van toetsing en monitoring;
-de onzekerheid in het Nederlands gemiddelde van de stikstofdepositie
bedraagt 20-30% (2𝜎𝜎).
De onzekerheidsschatting op lokale schaal is qua ruimtelijke resolutie
representatief voor een willekeurig meetlocatie en de omgeving direct
rondom dat meetpunt: de zogeheten ‘footprint’. Dit is het gebied
bovenwinds dat invloed heeft op de meting, ongeveer 100-1000 m. De
ruimtelijke resolutie van DN (hexagoon van 1 ha) als de GDN (grid-cel
van 1 × 1 km2) liggen dus beide in de range van deze
representativiteitsschatting. Op basis van dit onderzoek kunnen we niet
onderscheiden of deze onzekerheid representatief is voor een ha of een
km2. Maar de kleinere schaal van DN heeft wel als voordeel dat variaties
in de depositie binnen 1 × 1 km2 zichtbaar worden. Dit kan, omdat
verschillen in de locatie van bronnen of lokale eigenschappen van het
aardoppervlak per hexagoon een andere depositiewaarde kunnen
opleveren.
Het effect van de resolutie op de onzekerheidsschatting kent een
tegenstelling. Aan de ene kant laat de vergelijking tussen DN en de GDN
zien dat er binnen een 1 × 1 km2 grid-cel forse verschillen kunnen
optreden. Met name als een bron en het meetpunt dicht bij elkaar
zitten. Dat wijst erop dat berekeningen op de fijne schaal zinvol zijn,
omdat dan variaties in de depositie tussen hexagonen zichtbaar zijn.
Anderzijds wordt een sterk dempend effect op de onzekerheid
waargenomen als de depositie over een groter gebied wordt bepaald.
Om dat te kwantificeren, is er meer en andersoortig onderzoek nodig.
Bij toetsing over grotere gebieden, zoals de Commissie Hordijk adviseert
(Hordijk et al., 2020), zou het dempend effect inderdaad tot kleinere
onzekerheden in de depositie kunnen leiden. Maar merk op dat dit dus
alleen geldt voor de gemiddelde depositie over het gebied.

Page 22 of 93

RIVM report 2022-0085

In dit statusrapport hebben we praktische keuzes over de uitvoering
moeten maken:
• In de zes depositiecomponenten is een aantal componenten
samengenomen. Zo zitten NOy, HNO3- en nitraathoudende
deeltjes in de droge depositie van secundaire componenten.
• In de schatting van de correlatie tussen de diverse onzekerheden
hebben we de feedback tussen depositiesnelheid en concentratie
niet meegenomen. Bij een analyse die daar meer op is
toegespitst, zou dit tot lagere onzekerheidsschattingen kunnen
leiden.
• De aanwezige informatie over depositiesnelheden van deeltjes en
stikstofoxiden is erg beperkt.
De diverse verbeterpunten komen aan bod in een volgend statusrapport.
Ook zullen we dan kijken naar de onzekerheden in de sectorbijdragen en
naar de onzekerheid van een individuele bron of projectbijdrage.
De systematische afwijkingen tussen modelberekeningen en metingen
worden gecorrigeerd door de modelberekeningen te kalibreren aan de
metingen. Daarmee verbetert de (nauwkeurigheid van de)
depositieberekening. Maar op locaties waar een sterke kalibratieslag
nodig is, zullen interpretaties van modeluitkomsten minder nauwkeurig
zijn.
Er zijn diverse onderzoeken geprogrammeerd binnen het NKS om de
inzichten in processen en onzekerheden te verbeteren en deze
systematische onzekerheden te verkleinen. Bijvoorbeeld door het doen
van meer metingen, of door vergelijkingen met andere modellen en een
gevoeligheidsanalyse van de diverse parameters in het OPS-model. Dit
sluit aan bij de ambitie van het ministerie van LNV om een transparante
en wetenschappelijk solide basis te hebben voor het nationale
stikstofbeleid met gekwantificeerde onzekerheden.

Page 23 of 93

RIVM report 2022-0085

Page 24 of 93

RIVM report 2022-0085

2

Introduction

2.1

Aim of the study
The deposition of nitrogen compounds may result in a loss of
biodiversity in nature areas. In the Netherlands, the deposition of
nitrogen is exceeding the critical loads on a large scale. 2 Therefore, the
current policy in the Netherlands focuses on reducing the nitrogen
deposition.

Figure 1 Reactive nitrogen in the air, left: the sources of reduced (NHₓ, black)
and oxidised (NOₓ, grey), right: dry (yellow) and wet (blue) deposition.

The nitrogen deposition consists of various compounds: oxidised
nitrogen (such as nitrogen oxides and nitrate) and reduced nitrogen
(such as ammonia and ammonium). The oxidised nitrogen originates
mainly from emissions from traffic and industrial activities by way of
combustion processes. The reduced nitrogen largely originates from
agricultural activities.
The nitrogen deposition in the Netherlands is determined by means of a
combination of calculations based on the OPS model (Sauter et al.,
2023; van Jaarsveld, 2004) and measurements. The OPS model is the
basis for two sets of calculations. First, RIVM calculates the total
(annual) nationwide nitrogen deposition at a 1 × 1 km2 resolution for
the GCN/GDN maps (Grootschalige Concentratiekaart Nederland/
Grootschalige Depositiekaart Nederland, Hoogerbrugge et al., 2022).
Second, RIVM calculates at a hectare scale for a subset of Natura 2000
areas for the DN maps (Depositiekaarten stikstof Natura 2000-gebieden,
https://monitor.aerius.nl/).
Uncertainty plays an important role in the debate on nitrogen deposition
numbers and their policy implications. Uncertainties are present in many
aspects of the current method by which the nitrogen deposition is
2
Critical load is defined as “the highest deposition of nitrogen as NOX and/or NHY below which harmful effects in
ecosystem structure and function do not occur according to present knowledge” (ICP, M&M, 2004).

Page 25 of 93

RIVM report 2022-0085

determined: in the emissions, in the modelling of the atmospheric
processes (dispersion, transport, chemistry, deposition) and in
measurements. To identify areas for improvement, frequent scientific
reviews have been carried out (Hordijk et al., 2020; Sutton et al.,
2015). In these reviews the importance of uncertainty analysis is
identified.
Analysis of the uncertainty in concentrations and wet deposition are
available in, for example, the EMEP status report (EMEP, 2022). The
EMEP does not report the uncertainty of the total nitrogen deposition.
Walker et al. (2019) state that complete uncertainty analysis of the
modelled total nitrogen deposition is not yet available.
The government aims at a transparent and scientifically robust basis for
the national nitrogen policy that quantifies uncertainties. To achieve
this, the Dutch Ministry of Agriculture, Nature and Food Quality
(ministerie van Landbouw, Natuur en Voedselkwaliteit: LNV) has
initiated a National Knowledge Programme on Nitrogen (Nationaal
Kennisprogramma Stikstof: NKS). The NKS programme will propose and
implement improvements in the measurement and modelling system
that is used to evaluate the nitrogen deposition on nature areas. Thus, a
thorough quantification of the uncertainties in determining the nitrogen
deposition is important. This is not only important for policy support and
communication purposes, but also for setting up a reference to quantify
the uncertainty reduction by improving the model and measurement
system. Moreover, the knowledge on the most important sources of
uncertainties will be used to prioritise research activities.
Insight into and quantification of uncertainties is wanted in various types
of (policy) indicators or calculations, including:
• uncertainty in the nitrogen deposition on the Netherlands and its
nature areas;
• uncertainty in the nitrogen deposition on 1 × 1 km2 and ha
scales;
• uncertainty in the exceedance of the critical loads of nature areas
and the area(s) where critical loads are exceeded;
• uncertainty in the trends in the nitrogen deposition in time;
• uncertainty in the emission contribution by sectors to the
nitrogen deposition on the Netherlands;
• uncertainty in the contribution by a single source to the nitrogen
deposition on a specific nature area on a local scale, i.e. within
25 km from the source;
• uncertainty in the contribution by a single source to the nitrogen
deposition on a specific nature area on a large scale, i.e. up to a
few hundred km.
In general, no detailed uncertainty analyses of the above indicators or
calculations have been carried out yet. However, Van Jaarsveld (2004)
did analyse the uncertainty in the deposition. On the basis of this
analysis, a 70% uncertainty was estimated in the first publication of the
GDN (Velders et al., 2010). A more precise estimation of the uncertainty
on the regional and national scales is desirable, as is an uncertainty
estimate on a local scale. In this report, we will estimate the
uncertainties in the nitrogen deposition on national and local scales.

Page 26 of 93

RIVM report 2022-0085

This report should be considered as a status report in which we give an
outline the current knowledge on the various uncertainties and options
to improve the uncertainty estimates.
On the basis of these options, along with new methods and
measurements becoming available from the NKS programme and
literature, we will update this report. Subsequently, we will also take
uncertainty estimates for other indicators or calculations into account.
2.2

Detailed background
In this section, we will give a more detailed overview of the various
N compounds, and of the processes and methods that are important in
this study.
As stated before, nitrogen deposition consists of both oxidised and
reduced nitrogen. Both groups contain primary compounds as well as
compounds that are formed by chemical reactions in the atmosphere
(secondary compounds). For instance, particles called secondary
inorganic aerosols (SIAs).
Oxidised nitrogen contains one or more oxygen atoms. The main species
in the gas phase are nitrogen oxides (NO, NO2 and HNO3), while in SIA,
it is nitrate (NO3‒).
Reduced nitrogen contains hydrogen atoms. The main species are
ammonia (NH3) in the gas phase and ammonium NH4+ as the secondary
compound.
The amount of these compounds per unit volume of air is called the
concentration. The amount of these compounds that get to the surface
per surface area (e.g. m2 or ha) per unit of time is called deposition. The
concentration influences the deposition. In general, the larger the
concentration, the larger the deposition on the Earth’s surface.
Two mechanisms of deposition can be distinguished (see Figure 1, righthand side). First, the compounds get to the surface in and with
precipitation, which is called wet deposition. Second, plants and soils
absorb nitrogen compounds from the air, which is called dry deposition.
Apart from concentration, there are other processes and factors that
influence deposition. These differ across the six compounds. For
instance, wet deposition is influenced by the amount and intensity of the
precipitation.
For dry deposition, the deposition velocity is important. The deposition
velocity is a measure of how easily the compounds move through the air
towards the Earth’s surface and are subsequently absorbed by soil and
plants. The deposition velocity differs for each N compound. And it also
strongly depends on the weather, land use and roughness of the Earth’s
surface, soil type and vegetation type.
Since the deposition velocity is important to the calculation of the dry
deposition, it is also important to have an estimate of uncertainty within
the deposition velocity. Therefore, it is important to note that while the
concentration influences the deposition, it also works the other way
around: the amount of deposition has an effect on the concentration
because mass is removed from the atmosphere by that deposition and
thus reduces the concentration. We will analyse this effect in Chapter 5.
Given the different species and mechanisms, in this report, we
Page 27 of 93

RIVM report 2022-0085

Page 28 of 93

RIVM report 2022-0085

distinguish six compounds that contribute to the total deposition:
𝐷𝐷N tot = 𝐷𝐷NH4 wet + 𝐷𝐷NH3 primary dry + 𝐷𝐷NH4 secondary dry +
𝐷𝐷NO3 wet + 𝐷𝐷NOx primary dry + 𝐷𝐷NO3 secondary dry

Eq. 1

For each compound, the model calculations are calibrated against
measurements. In this way, the best estimate on the levels of nitrogen
deposition on the Netherlands is obtained. For the calculations of the
wet deposition, we calibrate the model against the observations of the
wet deposition. In the case of dry deposition, insufficient measurement
data is available for direct calibration. Moreover, there are no real-time
measurements of the deposition velocity. Therefore, the relation
between the modelled and measured air concentration of the
compounds is used to calibrate the dry deposition.
In Chapter 3 and 4, we will calculate the uncertainty of both the
uncalibrated and the calibrated values per compound and the overall
uncertainty.
2.3

Reading Guide
In Chapter 3, we will present the results on the uncertainty estimates in
the various deposition compounds. Subsequently, in Chapter 4, we will
describe the overall uncertainty in the total nitrogen deposition. In
Chapter 5, we will present the effect on uncertainty in the nitrogen
deposition due 6to the feedback mechanism of the deposition processes
on the concentration. Then in Chapter 6 we will show the comparison
between the deposition maps on 1 x 1 km2 (GCN/GDN) and ha scales
(DN). We will finish off with the conclusion and recommendations in
Chapter 7.

Page 29 of 93

RIVM report 2022-0085

Page 30 of 93

RIVM report 2022-0085

3

Uncertainty in nitrogen deposition estimates per compound

Authors: Koen Siteur, Joost Wesseling, Ronald Hoogerbrugge, Cor
Jacobs
In this chapter, we estimate the uncertainties of the six compounds of
the nitrogen deposition by comparing OPS model outputs to
measurements over the 2005-2021 period. We evaluate the residual
uncertainty. So we define the uncertainty as the residuals after
calibration of the model output against measurements and that are
assumed to be random.
3.1

Measurements
The aim of this chapter is to estimate the uncertainty of the modelled
deposition fluxes by comparing the model results to measurements.
However, deposition measurements are not available for each and every
compound. This means that we partly rely on concentration
measurements. The measurements used in this study are listed in
Table 1. All the measurements are annual averages.
The dataset consists of measurements used in the production of
GCN/GDN maps, combined with the data collected by the relatively
recent and expanding network of COTAG systems that measure the dry
deposition of NH3 (Rutledge-Jonker et al., 2023).
Table 1 also lists the uncertainty estimates for each of the measurement
methods, which were derived in previous studies. The listed estimates
relate to the uncertainty of the annual average of the measured
quantity. NO2 measurements by passive samplers are performed in
duplicates to reduce measurement uncertainty. For NH3, at least one
site per nature area has triplicate measurements. To account for
duplicates and triplicates, the measurement uncertainty of the passive
samplers of NO2 and NH3 is divided by the square root of the number of
passive samplers per site.
For most of the compounds, an approximately equal number of
measurements is available for each of the years included in this study.
Two exceptions are NH3 and NO2 concentrations, for which the number
of sites has been increasing in recent years. For NO2, the number of
measurement sites has doubled in 2019, to approximately ninety sites,
by including measurements of Palmes diffusion tubes performed in the
‘Samen Meten’ programme and the MAN network (Siteur et al., 2021).
For NH3, the number of sites in the MAN network has been increasing
steadily since 2005 (man.rivm.nl). As a result, the uncertainty estimates
that are derived in this chapter are slightly biased towards recent years.

Page 31 of 93

RIVM report 2022-0085

Table 1 Measurements included in this study. 𝑛𝑛 is total number of measurements included in this study, and the number of
measurements performed in 2021 has been added between brackets. The uncertainty column provides an equation of the
uncertainty of the measurement (𝑠𝑠 2 ) because, in most cases, the uncertainty depends on the flux (𝐹𝐹) or on the concentration (𝐶𝐶).

Quantity
Dry deposition
NH3
Concentration NH3

Unit
mol/ha/y

Method/device
COTAG

Period
2012-2020

n
18 (3)

μg/m3

Gradko tubes

2005-2021

Concentration
NH4+

μg/m3

2005-2008*
2009-2012
2013-2021**

Wet deposition
NHx
Concentration NO2

mol/ha/y

Low Volume
Sampler
Leckel
Leckel
Wet only samplers

3123
(278)
18
16
36 (4)

2005-2021

135 (6)

μg/m3

Chemiluminescence
Palmes tubes

2005-2021
2017-2021

651 (41)
142 (51)

Concentration
NO3-

μg/m3

2005-2008*
2009-2012
2013-2021**

18
16
36 (4)

Wet deposition
NOy

mol/ha/y

Low Volume
Sampler
Leckel
Leckel
Wet only samplers

2005-2021

135 (6)

Uncertainty
𝑠𝑠 2 = (0.1𝐹𝐹)2

Reference
(Rutledge-Jonker et al.,
2023)
(Noordijk et al., 2020)

𝑠𝑠 2 = (0.080𝐶𝐶)2
𝑠𝑠 2 = (0.080𝐶𝐶)2
𝑠𝑠 2 = (0.101𝐶𝐶)2

Stefess, 2022, pers.
com.

𝑠𝑠 2 = 0.302 + (0.06𝐶𝐶)2

𝑠𝑠 2 = (0.031𝐹𝐹)2

𝑠𝑠 2 = 0.9932 + (0.0412𝐶𝐶)2
𝑠𝑠 2 = (0.121𝐶𝐶)2
𝑠𝑠 2 = (0.052𝐶𝐶)2
𝑠𝑠 2 = (0.052𝐶𝐶)2
𝑠𝑠 2 = (0.070𝐶𝐶)2

𝑠𝑠 2 = (0.031𝐹𝐹)2

Stefess, 2022, pers.
com.
Stefess, 2022, pers.
com.
(Nguyen and Wesseling,
2016)
Stefess, 2022. pers.
com.
Stefess, 2022, pers.
com.

*For the Low Volume Sampler used until 2008, no uncertainty estimates are available. The used uncertainty estimate is probably an underestimation.
**For the 2013-2021 period, NH4+ and NO3- concentrations were measured every other day, increasing the uncertainty in the annual average
concentrations.

Intern gebruik

RIVM report 2022-0085

3.2

Model calculations
We performed model runs of OPS-LT version 5.1.1.0 to obtain modelled
annual average concentrations and deposition fluxes at the
measurement sites over the 2005-2021 period. These runs on receptor
points provide the concentrations and deposition fluxes at the exact
coordinates of the measurement sites. The model input consists of
annual NH3 and NOx emissions at a 1 x 1 km2 resolution, with a spatial
distribution of each GCN/GDN emission sector set to the 2020 value.
The sums of the Dutch emissions were scaled per sector and per year,
following the emission data of the Emission Registration (ER;
www.emissieregistratie.nl). Emission sums of foreign sources were
scaled using emission data of the Centre on Emission Inventories and
Projections (CEIP; www.ceip.at). The meteorological statistics that we
used to run OPS were calculated for each year using METPRO with
meteo data from the Royal Netherlands Meteorological Institute (KNMI;
www.knmi.nl) as input. Land use and roughness maps were derived
from LGN2020 (www.lgn.nl).
The resolution of the land use and roughness maps determines the
resolution of the dry deposition calculations. As GDN maps are currently
produced at a resolution of 1 x 1 km2, the analysis in this chapter is
based on 1 x 1 km2-resolution land use and roughness maps. We realise
that this resolution may not match the footprint of the COTAG dry
deposition measurements, making a comparison with measurements
less reliable for estimating the uncertainty. To study the sensitivity of
this choice, we also perform model runs at resolutions of 0.25 x 0.25,
0.5 x 0.5, 2.5 x 2.5 and 5 x 5 km2.

3.2.1

Post-processing of OPS outputs
OPS produces outputs of concentrations and deposition fluxes for each
compound. OPS also computes annual average dry deposition velocities,
which, however, are not included in the model outputs. As the dry
deposition velocity is required for uncertainty estimation of the dry
deposition fluxes (see Section 3.5.1), a post-processing step is required.
The dry deposition velocity is obtained by dividing the annual, modelled
deposition flux by the modelled annual concentration:
𝑉𝑉d,X =

𝐹𝐹𝑑𝑑,𝑋𝑋
𝐶𝐶𝑋𝑋

Eq. 2

where 𝐶𝐶 is the concentration of compound 𝑋𝑋, 𝐹𝐹d is the dry deposition flux
of that compound and 𝑉𝑉d is the dry deposition velocity. Note that OPS
takes into account that the vegetation can (re-)emit NH3 to the
atmosphere, by using a compensation point for the vegetation. Hence,
for NH3, Eq. 2 yields the effective annual averaged deposition velocity.
OPS calculates and returns the concentrations of all the measured
compounds listed in Table 1, except for NO2. Within OPS, the
NO2 concentrations used to compute the NO2 deposition flux are derived
from the NOx concentration with the following empirical relationship
(Sauter et al., 2023):
𝐶𝐶NO2 = 𝛽𝛽1 ln�𝐶𝐶NOx � + 𝛽𝛽2
𝛽𝛽

𝛽𝛽

for 𝐶𝐶NOx > exp �1 − 𝛽𝛽2 �
1

𝛽𝛽

𝐶𝐶NO2 = 𝛽𝛽1 exp �𝛽𝛽2 − 1� 𝐶𝐶NOx for 𝐶𝐶NOx < exp �1 − 𝛽𝛽2 �
1

Eq. 3

1

Page 33 of 93

RIVM report 2022-0085

where 𝛽𝛽1 = 8.6 and 𝛽𝛽2 = −12.4 (concentrations in ppb). As this
relationship is used within OPS, we also use it in the computation of the
NO2 concentration from the NOx concentration.
Finally, in the model output, the dry deposition fluxes are grouped into
the deposition of primary and secondary compounds.
For the reduced compounds, the primary compound is NH3, and the
secondary compound is NH4+.
For the oxidised compounds in the model, the primary compounds are
NOx = NO + NO2 + HNO2, and the secondary compounds are NO3- +
HNO3.
In this study, we derive relative uncertainties in concentrations and
ultimately deposition fluxes of NO2 and NO3, and we assume that these
are representative of the uncertainties of the primary and secondary
compounds, respectively.
3.3

Removing systematic model bias
The biases of the model are corrected by calibrating the modelled
deposition fluxes against measurements. We use the calibration
techniques as applied in constructing the GDN maps. That is, the wet
deposition fluxes are calibrated against measured deposition fluxes
(Hoogerbrugge et al., 2022) and the dry deposition of NH3 is calibrated
against NH3 concentration measurements of the MAN network (Wichink
Kruit et al., 2020). In Table 2, these calibration variables are marked in
green.
The dry deposition fluxes of NH4+, NO3-, and NO2, are currently not
calibrated in the GDN procedure. To estimate the bias in the deposition
fluxes of these compounds, the modelled concentrations are calibrated
against measured concentrations using the methods currently applied to
calculate the official GCN concentration maps. In Table 2, these
calibration variables are marked in blue.
Biases in the dry deposition velocities of NH3, NH4+, NO2 and NO3- could
not be estimated. In Table 2, these dry deposition velocities are marked
in orange.
Table 2 Calibration methods to estimate biases in the deposition fluxes. Column
3 and 4 indicates if the variables are calibrated: N means No and Y means Yes.
Column 3 (T) for the analysis in This study and 4 (G) in the GDN procedure.

Flux

𝐹𝐹d,NH3 = 𝑉𝑉d,NH3 ⋅ 𝐶𝐶NH3

𝐹𝐹d,NOx =
𝑉𝑉d,NOx
⋅ �𝐶𝐶NO + 𝐶𝐶NO2 + 𝐶𝐶HNO2 �

𝐹𝐹d,NH+4 = 𝑉𝑉𝑑𝑑,𝑁𝑁𝑁𝑁4+ ⋅ 𝐶𝐶𝑁𝑁𝐻𝐻4
𝐹𝐹d,HNO +NO− =
3

3

𝑉𝑉d,HNO +NO− ⋅ (𝐶𝐶NO3
3

3

+ 𝐶𝐶HNO3 )
𝐹𝐹w,NHx
𝐹𝐹w,NOy

Page 34 of 93

𝑉𝑉

Variable
𝑉𝑉d,NH3
𝐶𝐶NH3
𝑉𝑉d,NOx
𝐶𝐶NO
𝐶𝐶NO2
𝐶𝐶HNO2
𝑉𝑉d,NH+4
𝐶𝐶NH4

d,HNO3 +NO−
3

𝐶𝐶NO3
𝐶𝐶HNO3

T
N
Y
N
N
Y
N
N
Y

G
N
Y
N
N
N
N
N
N

Method/Reference

N
Y
N

N
N
N

Linear regression through origin

Y
Y

Y
Y

Linear regression through origin
Linear regression through origin

Wichink Kruit et al., 2020

Siteur et al., 2021

Linear regression through origin

RIVM report 2022-0085

3.4

Uncertainty in the model
The combined uncertainty in the model and measurements is obtained
by:
2
𝑠𝑠OPS
=

2
∑𝑛𝑛
𝑖𝑖=1(𝑋𝑋𝑖𝑖,obs −𝑋𝑋𝑖𝑖,OPS )

𝑛𝑛

Eq. 4

where 𝑋𝑋𝑖𝑖,obs is one of the measured quantities in Table 1, and 𝑋𝑋𝑖𝑖,OPS is the
same quantity as modelled by OPS, both at measurement site 𝑖𝑖. The
number of measurements is represented by 𝑛𝑛.
Similarly, the uncertainty in the calibrated model output is obtained
using:
2
𝑠𝑠cal
=

∑𝑛𝑛𝑖𝑖=1(𝑋𝑋𝑖𝑖,obs − 𝑋𝑋𝑖𝑖,cal )2
𝑛𝑛

Eq. 5

where 𝑋𝑋𝑖𝑖,cal is the calibrated OPS output at measurement site 𝑖𝑖. In Eq. 5,
the modelled quantity 𝑋𝑋 at location 𝑖𝑖 is calibrated against all the
measurements except the measurement at location 𝑖𝑖. This is referred to
as leave-one-out cross-validation. This validation enables a fair
comparison between calibrated model and measurement, because the
measurements used for calibration are not used for validation (Siteur et
al., 2021).
The uncertainty obtained through Eq. 5 is composed of model and
measurement uncertainty. Random uncertainties in measurements have
been estimated in previous studies and are listed in Table 1. By
subtracting the random measurement uncertainty, 𝑠𝑠obs , from the
uncertainty estimated by means of Eq. 5, the uncertainty that can be
attributed to the model is obtained:

3.5

2
2
2
𝑠𝑠mod
= 𝑠𝑠cal
− 𝑠𝑠obs

Eq. 6

Uncertainties in dry deposition
Measurements of dry deposition fluxes are available for NH3 only, which
means that a direct uncertainty estimate of the modelled dry deposition
flux can only be made for NH3, using the procedure described in the
previous section. For the dry deposition fluxes of the other compounds,
no uncertainty estimate can be made directly. For these compounds, we
will approximate the uncertainty in the deposition flux by combining our
uncertainty estimate of the concentration 𝐶𝐶 with an approximation of the
uncertainty in the dry deposition velocity 𝑉𝑉𝑑𝑑 , as we will describe in this
section.

3.5.1

Decomposition of the uncertainty in dry deposition fluxes
Dry deposition fluxes can be calculated with the product of the
deposition velocity and the concentration:
𝐹𝐹d = 𝑉𝑉d ∙ 𝐶𝐶

Eq. 7

Page 35 of 93

RIVM report 2022-0085

Using a Taylor expansion, the uncertainty in 𝐹𝐹𝑑𝑑 can be approximated by:
𝜕𝜕𝐹𝐹𝑑𝑑 2 2
𝜕𝜕𝐹𝐹 2
𝜕𝜕𝐹𝐹 𝜕𝜕𝐹𝐹
� 𝑠𝑠𝑉𝑉𝑑𝑑 + � 𝑑𝑑 � 𝑠𝑠𝐶𝐶2 + 2 𝑑𝑑 𝑑𝑑 𝑠𝑠𝑉𝑉2𝑑𝑑𝐶𝐶 + ⋯
𝜕𝜕𝑉𝑉𝑑𝑑
𝜕𝜕𝜕𝜕
𝜕𝜕𝑉𝑉𝑑𝑑 𝜕𝜕𝜕𝜕
2
𝑠𝑠𝑉𝑉
𝑠𝑠𝑉𝑉
𝑠𝑠 2
𝑠𝑠
= 𝐹𝐹𝑑𝑑2 � 2𝑑𝑑 + 𝐶𝐶2 + 2 ⋅ 𝑑𝑑 ∙ 𝐶𝐶 ⋅ 𝑠𝑠𝑉𝑉2𝑑𝑑𝐶𝐶 + ⋯ �
𝑉𝑉
𝐶𝐶
𝑉𝑉
𝐶𝐶

𝑠𝑠𝐹𝐹2𝑑𝑑 = �

Eq. 8

𝑑𝑑

𝑑𝑑

where 𝑠𝑠𝑉𝑉2d 𝐶𝐶 = 𝜌𝜌(𝑉𝑉d , 𝐶𝐶)𝑠𝑠𝑉𝑉d 𝑠𝑠𝐶𝐶 , with 𝜌𝜌 being the correlation coefficient, 𝑠𝑠𝐶𝐶 is
the uncertainty in the concentration estimated using Eq. 5, and 𝑠𝑠𝑉𝑉d
denotes the uncertainty in the deposition velocity. The “...” denotes
higher-order terms.
In Appendix B, we will show that, on the local scales considered here,
the contribution to the covariance term is negligible. Assuming, 𝑉𝑉𝑑𝑑 and 𝐶𝐶
are indeed independent variables, the exact uncertainty in 𝐹𝐹d can be
written as (see Appendix B.1 and Goodman, 1960)
2
𝑠𝑠𝑉𝑉

𝑠𝑠 2

𝑠𝑠𝐹𝐹2𝑑𝑑 = 𝐹𝐹𝑑𝑑2 � 2𝑑𝑑 + 𝐶𝐶2 +
𝑉𝑉𝑑𝑑

3.5.2

𝐶𝐶

2
2
𝑠𝑠𝑉𝑉
𝑑𝑑 𝑠𝑠𝐶𝐶

𝑉𝑉𝑑𝑑2 𝐶𝐶 2

�

Eq. 9

which is the relationship that is used in the present analysis to estimate
the uncertainty in the dry deposition of NH4+, NO3-, and NO2. Note that
the uncertainty in 𝐹𝐹d is proportional to 𝐹𝐹d , and that the relative
uncertainty in 𝐹𝐹d is composed of the relative uncertainties in 𝑉𝑉d and 𝐶𝐶.

Estimation of uncertainty in dry deposition velocities of NH4+, NO3-, and
NO2
The uncertainty in 𝐹𝐹d can be computed for NH4+, NO3-, and NO2 using Eq.
9, with the relative uncertainty in C and an approximation of the relative
uncertainty in 𝑉𝑉d . For the present analysis, we will assume that the
relative uncertainty in 𝑉𝑉d can be approximated by the relative
uncertainty in 𝑉𝑉d for NH3. This is achieved by computing the uncertainty
in the effective dry deposition velocity of NH3 using dry deposition
measurements, and subtracting uncertainty resulting from uncertainties
in the estimation of the compensation point, as described below.
The dry deposition flux of NH3 is calculated with the product of the
effective deposition velocity and the NH3 concentration,
𝐹𝐹d = 𝑉𝑉d,eff ∙ 𝐶𝐶

Eq. 10

𝑉𝑉d,eff = 𝑓𝑓(𝐶𝐶)𝑉𝑉d

Eq. 11

where the effective deposition velocity accounts for deposition and
(re-)emission. The effective deposition velocity is determined by the
compensation concentration, which in turn is a function of the
atmospheric background concentration 𝐶𝐶:

Using these two relationships, we can derive a relative uncertainty in 𝑉𝑉d .
The uncertainty in the effective deposition velocity is calculated with:
2
𝑠𝑠𝑉𝑉

d,eff
2
𝑉𝑉d,eff

2
𝑠𝑠𝐹𝐹

𝑠𝑠 2

𝑠𝑠 2

= � 2d − 𝐶𝐶2� �1 + 𝐶𝐶2�

Page 36 of 93

𝐹𝐹d

𝐶𝐶

𝐶𝐶

−1

Eq. 12

RIVM report 2022-0085

This relationship can be obtained by rewriting the equation for exact
variance (Goodman, 1960). Here, 𝑠𝑠𝐹𝐹2𝑑𝑑 is estimated using COTAG
measurements, while 𝑠𝑠𝐶𝐶2 is estimated using the concentration
measurements from the MAN network.
Following the same reasoning, the uncertainty in the dry deposition
velocity, 𝑠𝑠𝑉𝑉2d can be obtained with:
2
𝑠𝑠𝑉𝑉

d
𝑉𝑉d2

2
𝑠𝑠𝑉𝑉

=�

d,eff
2
𝑉𝑉d,eff

−

2
𝑠𝑠𝑓𝑓(𝐶𝐶)

𝑓𝑓(𝐶𝐶)2

� �1 +

2
𝑠𝑠𝑓𝑓(𝐶𝐶)

𝑓𝑓(𝐶𝐶)2

−1

�

Eq. 13

This is the uncertainty that is adopted for the dry deposition velocities of
NH4+, NO3-, and NO2.
To get an estimate for the uncertainty resulting from the compensation
concentration, we need to define function 𝑓𝑓(𝐶𝐶). Both in the OPS model
and in observations over peatland (Cape et al., 2008), we find that for
relevant ranges of 𝐶𝐶, 𝑓𝑓(𝐶𝐶) can be approximated by:
𝑓𝑓(𝐶𝐶) = exp(−𝑎𝑎𝑎𝑎)

We estimate the exponent to be 𝑎𝑎 ≈ 0.16

Eq. 14

𝑚𝑚3
𝜇𝜇𝜇𝜇

, which is based on

observations by Cape et al. (2008) over peatland, across a
concentration range of 1 to 10 μg/m3.
Given the exponential function, the relative uncertainty in 𝑓𝑓(𝐶𝐶) becomes:
2
𝑠𝑠𝑓𝑓(𝐶𝐶)

𝑓𝑓(𝐶𝐶)2

= 𝑎𝑎2 𝑠𝑠𝐶𝐶2

Eq. 15

which can be substituted in Eq. 13 to obtain the uncertainty in the
deposition velocity. Note that 𝑠𝑠𝐶𝐶 is the absolute uncertainty in the
atmospheric background concentration. In OPS, this atmospheric
background concentration is obtained by iteratively running the model
until a constant concentration is reached. No spatial calibration is
performed in this iteration. Therefore, 𝑠𝑠𝐶𝐶 is the uncertainty in the
uncalibrated OPS result, which we obtain using the MAN NH3
concentration measurements.
3.6

Results
To illustrate the performance of the OPS model, Figure 2 and Figure 3
show the modelled and observed values of the measured quantities over
the whole observation period (2005-2021), both before and after
calibration. The calibration, which was performed for each separate year
using all observations except for the observation under consideration
(leave-one-out cross-validation), generally moves the modelled results
closer to the 1:1 line, without removing the scatter. This indicates that
the calibration removes much of the bias, without affecting the random
uncertainty. A notable exception is the NH3 concentration, for which the
calibration removes a significant amount of the scatter. This is because
part of the scatter in NH3 concentrations is caused by an east-west
gradient in the residuals, which is removed using a spatially explicit
calibration method (Wichink Kruit et al., 2020). The average bias in the
NH3 calculations is small.
Page 37 of 93

RIVM report 2022-0085

Figure 2 Modelled and observed concentrations and wet deposition fluxes of
oxidised compounds at the observation sites, before and after calibration.

Page 38 of 93

RIVM report 2022-0085

Figure 3 Modelled and observed concentrations and wet deposition fluxes of
reduced compounds at the observation sites, before and after calibration.

Page 39 of 93

RIVM report 2022-0085

Figure 4 Modelled and observed NH3 dry deposition fluxes (BAR: Bargerveen,
OVP: Oostelijke Vechtplassen, VEL: Hoge Veluwe). In the right-hand figure,
concentrations are calibrated, using the MAN observations, following the current
GCN/GDN procedure. Calibration of the resulting deposition flux is not yet
possible (see text).

Figure 4 shows measured and modelled dry deposition fluxes of NH3 at
the COTAG locations.
Table 3 summarises the model performance for each of the measured
quantities, both before and after calibration. The uncalibrated model
performs relatively well on NO2 concentrations and wet deposition
fluxes, as indicated by the low relative standard deviations. The
uncertainties in modelled NH3 concentrations and the concentrations of
the secondary compounds are relatively high. However, the calibration
improves the model performance on these concentrations considerably.
Table 3 Model performance before and after calibration.
Subscripts obs, OPS and cal refer to statistics of observations, uncalibrated
results and calibrated results, respectively. 𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 is the model uncertainty, i.e.
the uncertainty that remains after subtracting measurement uncertainty. The
operator 〈 〉 means the average over all results at the observation sites.

compound
NH3
NO2
NH4+
NO3NHx
NOy
NH3

quantity

unit

〈𝑋𝑋𝑜𝑜𝑜𝑜𝑜𝑜 〉

observed

uncalibrated

〈𝑋𝑋𝑐𝑐𝑐𝑐𝑐𝑐 〉

𝐶𝐶
𝐶𝐶
𝐶𝐶
𝐶𝐶
𝐹𝐹𝑤𝑤
𝐹𝐹𝑤𝑤
𝐹𝐹𝑑𝑑

µg/m3
µg/m3
µg/m3
µg/m3
mol/ha/y
mol/ha/y
mol/ha/y

5.2
18.1
1.7
4.6
391
194
604

5.7
19.0
2.0
3.9
373
215
428

5.2
18.1
1.7
4.5
384
193
403

Page 40 of 93

𝑠𝑠𝑜𝑜𝑜𝑜𝑜𝑜

〈𝑋𝑋𝑂𝑂𝑂𝑂𝑂𝑂 〉

𝑠𝑠𝑂𝑂𝑂𝑂𝑂𝑂

𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐

0.4
1.2
0.2
0.3
12
6
62

2.1
2.4
0.7
1.0
61
34
228

1.2
1.9
0.6
1.0
64
25
249

Calibrated

𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚
Calibrated
without 𝑠𝑠𝑜𝑜𝑜𝑜𝑜𝑜
1.1
1.4
0.6
1.0
62
24
241

RIVM report 2022-0085

Table 4 Relative uncertainties per deposition flux.

compound
NH3
NO2
NH4+
NO3NHx
NOy

quantity
𝐹𝐹𝑑𝑑
𝐹𝐹𝑑𝑑
𝐹𝐹𝑑𝑑
𝐹𝐹𝑑𝑑
𝐹𝐹𝑤𝑤
𝐹𝐹𝑤𝑤

𝑠𝑠𝑉𝑉𝑑𝑑
𝑠𝑠𝐹𝐹,𝑐𝑐𝑐𝑐𝑐𝑐
〈𝐹𝐹𝑐𝑐𝑐𝑐𝑐𝑐 〉
〈𝑉𝑉𝑑𝑑 〉
Calibrated
23% 56%* 62%
11%
46% 47%
36%
46% 60%
22%
46% 52%
- 17%
- 13%

𝑠𝑠𝐶𝐶,𝑐𝑐𝑐𝑐𝑐𝑐
〈𝐶𝐶𝑐𝑐𝑐𝑐𝑐𝑐 〉

𝑠𝑠𝐶𝐶,𝑚𝑚𝑚𝑚𝑚𝑚 𝑠𝑠𝑉𝑉𝑑𝑑,𝑚𝑚𝑚𝑚𝑚𝑚
𝑠𝑠𝐹𝐹,𝑚𝑚𝑚𝑚𝑚𝑚
〈𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚 〉 〈𝑉𝑉𝑑𝑑,𝑚𝑚𝑚𝑚𝑚𝑚 〉
〈𝐹𝐹𝑐𝑐𝑐𝑐𝑐𝑐 〉
Calibrated without 𝑠𝑠𝑜𝑜𝑜𝑜𝑜𝑜
21% 55%*
60%
8%
44%
45%
35%
44%
58%
21%
44%
50%
16%
13%

*For NH3, the uncertainties in the dry deposition velocity represent uncertainties in the
effective deposition velocity.

Relative uncertainties were computed by dividing the uncertainties listed
in Table 3 by the respective means. These relative uncertainties were
then used to calculate the relative uncertainties of the six deposition
fluxes using Eq. 9, as shown in Table 4.
Table 3 and Table 4 show that the relative uncertainties in the dry
deposition fluxes are large, compared to those in the wet deposition
fluxes. This can be attributed to the relatively large uncertainty in the
dry deposition velocities. Subtracting the measurement uncertainty has
little effect, which means that most of the uncertainty is caused by the
model, particularly the calculation of the deposition velocity.
3.7

Interpretation of difference between modelled and measured dry
deposition fluxes
The uncertainty estimates of the dry deposition fluxes are based on a
small number of observations at the three COTAG sites only. As the
difference between observations and the model is relatively large at
these sites (see Figure 4), our uncertainty estimates for the dry
deposition fluxes are relatively high. Figure 4 suggests that the model
mostly underestimates the dry deposition fluxes at the COTAG sites. To
investigate the cause of this apparent bias, we used NH3 concentrations
measured at the same sites (at ~5m above the surface) to compute the
effective dry deposition velocities. We compare these with the modelled
effective dry deposition velocities. From Figure 5, we conclude that the
difference between measured and modelled deposition fluxes can be
explained by an underestimation of the dry deposition velocities at the
COTAG sites. As we only use this limited dataset, it is impossible to
assess if this bias is a systematic model bias that also applies to other
areas in the Netherlands. Therefore, we have compared the modelled
and measured deposition velocities with measurements reported in the
literature and with observations with a high resolution deposition flux
equipment (GRAHAM) at Solleveld (Vendel et al., 2023), as described in
Appendix E.1. On the basis of these comparisons, we cannot confirm
that OPS systematically underestimates the dry deposition velocity.
Therefore, no bias correction is applied to the modelled dry deposition
velocities in this study. A more detailed analysis of the differences
between the calculated and measured dry deposition fluxes will be
carried out. The improvement on the knowledge on deposition velocities
is an important objective in the Dutch research programme ‘Nationaal
Kennisprogramma Stikstof’ (NKS).
Page 41 of 93

RIVM report 2022-0085

Figure 5 Modelled and measured NH3 effective dry deposition velocities (left),
and concentrations (right) at the COTAG sites (BAR: Bargerveen, OVP:
Oostelijke Vechtplassen, VEL: Hoge Veluwe). The modelled deposition velocities
are obtained by dividing modelled deposition fluxes by modelled concentrations
at a height of ~5m. The observed deposition fluxes are obtained by dividing
measured deposition fluxes (COTAG) with measured concentrations (MAN) at a
height of ~5m.

3.8

Conclusion & Discussion
In this chapter, we have derived estimates of the uncertainties in six
modelled deposition fluxes. The dry deposition fluxes have the largest
relative uncertainty, mainly due to uncertainties in the dry deposition
velocities. Both modelled and observed deposition velocities for
ammonia fall well within the range of values reported in the
international literature. Measurement uncertainty contributes little to the
total uncertainty, implying that the differences between modelled and
measured fluxes and concentrations are mainly caused by model
uncertainties, particularly uncertainties in the deposition velocities.
Given the small number of sites for which dry deposition measurements
are available, it is hard to get reliable uncertainty estimates of the dry
deposition fluxes. This can be solved in future studies by expanding the
number of COTAG sites, and/or by deriving uncertainty estimates of the
dry deposition velocities using uncertainties in underlying variables, such
as the friction velocities, land use and roughness lengths. Finally, for a
fair comparison between model and measurement, the resolution of the
roughness and land use maps in the model should match the footprint of
the COTAGs. In present analysis, we have opted for a resolution of
1 × 1 km2, as this is the resolution currently used in the production of
GDN maps. A more detailed analysis of the differences between the
calculated and measured dry deposition fluxes at the COTAG-sites needs
to be made.

Page 42 of 93

RIVM report 2022-0085

4

Uncertainty in total nitrogen deposition

Authors: Koen Siteur, Joost Wesseling, Ronald Hoogerbrugge
In Chapter 3, we derived the relative uncertainties in each of the six
modelled deposition fluxes. However, the uncertainty of the six
individual fluxes does not yet inform us about the uncertainty in the
total nitrogen deposition. In addition, the derived uncertainties
represent the performance of calibrated model fluxes. In GDN and DN
maps, not all fluxes are currently calibrated, and the bias should be
accounted for in our uncertainty estimate. In addition, the uncertainty of
the uncalibrated model result is of interest, as it informs us about the
performance of the OPS model.
In this chapter, we establish a way to compute the uncertainty in the
total nitrogen deposition, given a combination of the six deposition
fluxes. We then compute the random uncertainty and bias for each
GCN/GDN grid cell. Finally, we determine the average uncertainty for
the total nitrogen deposition on Natura 2000 habitats and approximate
the uncertainty in the national average total nitrogen deposition.
4.1

Uncertainty in total nitrogen deposition
Since the six deposition fluxes are linked through chemical and physical
processes, both in the model and in the real world, they cannot be
viewed as independent random variables.
The variance in the sum of dependent random variables is calculated
with the sum of (co-)variances. Hence, the uncertainty in the total
nitrogen deposition is:
6

5

6

𝑖𝑖=1

𝑗𝑗=1 𝑖𝑖=𝑗𝑗+1

2
𝑠𝑠𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹
= � Var(𝐹𝐹𝑖𝑖 ) + 2 � � Cov�𝐹𝐹𝑖𝑖 , 𝐹𝐹𝑗𝑗 �

Eq. 16

The first term in Eq. 16 is the sum of the squared random uncertainties
2
, as derived in Chapter 3. The
of each of the six deposition fluxes, ∑6𝑖𝑖=1 𝑠𝑠𝐹𝐹𝐹𝐹
second term results from covariance between uncertainties. The
covariances in the second term can be written as:
Cov�𝐹𝐹𝑖𝑖 , 𝐹𝐹𝑗𝑗 � = 𝑠𝑠𝐹𝐹𝑖𝑖 𝑠𝑠𝐹𝐹𝑗𝑗 𝜌𝜌(𝐹𝐹𝑖𝑖 , 𝐹𝐹𝑗𝑗 )

Eq. 17

Where 𝜌𝜌 is the correlation coefficient between flux pairs. In Appendix
C.2, we use a result from Bohrnstedt and Goldberger (1969) to establish
that 𝜌𝜌 is solely determined by the relative uncertainties in 𝐶𝐶, 𝑉𝑉𝑑𝑑 and 𝐹𝐹𝑑𝑑 ,
as derived in Chapter 3, and by the correlation between the
concentrations and the dry deposition velocities. Relative uncertainties
being constant, this implies that 𝜌𝜌 is constant, i.e. it is independent of
the model result to which we wish to apply Eq. 16. In Appendix C.3, we
will compute the correlation matrices for 𝑉𝑉𝑑𝑑 and 𝐶𝐶 using residuals at the
measurement sites, supplemented with model results. These correlation
matrices are combined to obtain the correlation matrix for 𝐹𝐹, as
presented in Table 5.
Page 43 of 93

RIVM report 2022-0085

Table 5 Correlation matrix for the six deposition fluxes. This matrix was derived
using the model uncertainties (𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 ) from Chapter 3. The correlation matrix
based on 𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐 is very similar (differences are of a magnitude of ≈0.01) and is
shown in Appendix C.3.

𝝆𝝆

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑

1.00

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−
𝟑𝟑

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

𝟑𝟑

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱

0.74
1.00

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

0.58
0.55
1.00

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

1.00

0.30
1.00

0.22
0.67
1.00

𝟑𝟑

0.66
0.57
0.76

𝟑𝟑

-0.03
-0.01
0.34

-0.10
-0.05
0.24

Next, we can compute the absolute uncertainty for each deposition flux
𝑠𝑠𝐹𝐹𝑖𝑖 for each GDN grid cell for each year. We do this by multiplying the
relative uncertainties derived in Chapter 3 by the (calibrated) modelled
fluxes. This, in combination with Table 5, is then used to calculate the
uncertainty in the total deposition flux, averaged over all years. This
procedure yields Figure 6, a map in which the uncertainty in total
nitrogen deposition is averaged over the 2005-2021 period.
Note that the analysis described above assumes that uncertainties in 𝐶𝐶
and 𝑉𝑉𝑑𝑑 are proportional to the modelled 𝐶𝐶 and 𝑉𝑉𝑑𝑑 . This assumption will
be evaluated in Appendix A by comparing the assumed proportional
variance model with binned absolute uncertainty estimates. From this
comparison we conclude that for the current analysis, the proportional
model can be applied.

Figure 6 Random model uncertainty (𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 ) in total nitrogen deposition, averaged
over the 2005-2021 period.

Page 44 of 93

RIVM report 2022-0085

4.2

Total uncertainty in total nitrogen deposition
The map derived in the previous section reflects the uncertainty after
removing the known model biases from the model result. This remaining
uncertainty we interpret as random uncertainty, and it contains both the
unbiased uncertainties in the model and uncertainties resulting from the
bias correction. To obtain an estimate of the total model uncertainty, we
again add the known model biases to the obtained uncertainty
estimates.

4.2.1

Sum of model biases
The bias for each flux 𝑖𝑖 is determined by subtracting the calibrated
model result from the uncalibrated model result for each grid cell:
Δ𝑖𝑖 = 𝐹𝐹OPS,𝑖𝑖 − 𝐹𝐹cal,𝑖𝑖

Eq. 18

The sum of the biases becomes the total model bias:
6

Δ = � Δ𝑖𝑖

Eq. 19

𝑖𝑖=1

Figure 7 shows the average total model bias over the 2005-2021 period,
for the OPS model (Δ𝑂𝑂𝑂𝑂𝑂𝑂 ) and for the GDN maps (Δ𝐺𝐺𝐺𝐺𝐺𝐺 ). The calibration
of the GDN maps, and specifically the NH3 concentration map, reduces
the known bias by a factor of ~50. However, the known bias in the GDN
maps is not zero, as not all measured variables are currently used in the
calibration of GDN maps (see Chapter 3, Table 2). The major spatial
signature of the estimated bias in the OPS model is strongly reduced by
the calibration.

Figure 7 Known bias in total nitrogen deposition of the OPS model and of the
partially calibrated GDN map. These maps are obtained by subtracting the
calibrated maps from uncalibrated maps (Eq. 18). Note the factor-ten difference
in the range of the colour bar.

The estimation of the bias is itself subject to uncertainty. Bias
corrections that are based on a large number of observations yield a low
uncertainty in the estimation of the bias, and vice versa. This
uncertainty is accounted for by cross-validation, which means that it is
Page 45 of 93

RIVM report 2022-0085

part of the uncertainties derived in Chapter 3. As we currently are
unable to assess a bias in the deposition velocities, Figure 7 should be
interpreted with care. The distribution of known biases suggests an
underestimation in the west of the Netherlands and an overestimation in
the east, implying that model processes or model inputs are less well
captured in these regions.
4.2.2

Estimation of total uncertainty in total nitrogen deposition
Now that we have established the random uncertainty and bias for each
GDN grid cell, we will estimate the total uncertainty by assuming that
the local random uncertainty and the local uncertainty due to
uncorrected bias are mutually independent. The total uncertainty of the
uncalibrated OPS output, can be estimated as:
𝑠𝑠𝐹𝐹𝑡𝑡𝑡𝑡𝑡𝑡,OPS = �𝑠𝑠𝐹𝐹2𝑡𝑡𝑡𝑡𝑡𝑡,𝑂𝑂𝑂𝑂𝑂𝑂 + Δ2OPS

Eq. 20

𝑠𝑠𝐹𝐹𝑡𝑡𝑡𝑡𝑡𝑡,GDN = �𝑠𝑠𝐹𝐹2𝑡𝑡𝑡𝑡𝑡𝑡,𝐺𝐺𝐺𝐺𝐺𝐺 + Δ2GDN

Eq. 21

The total uncertainty in the GDN maps is becomes:

Figure 8 presents the resulting maps. Note that, since the estimated
known bias is smaller than the estimated random uncertainty, it
contributes little to the total uncertainty.

Figure 8 Sum of the random model uncertainty and the bias (𝛥𝛥) for OPS and
GDN.

4.3

Average uncertainty in nitrogen deposition on Natura 2000
habitats
The maps in the previous sections show estimates of the uncertainties
for each 1 × 1 km2 GDN grid cell. Since we are interested in the
uncertainty of the total nitrogen deposition on sensitive Natura 2000
habitats, we computed the average uncertainty, weighted by the
fractional cover of sensitive Natura 2000 habitats (Figure 9). A
Natura 2000 habitat is considered sensitive if its critical nitrogen

Page 46 of 93

RIVM report 2022-0085

deposition load is less than 2400 mol/ha/y (Dobben and Hinsberg,
2008).

Figure 9 Fractional cover of relevant Natura 2000 habitat for all 1 × 1 km2 GDN
grid cells.

Table 6 shows the average model uncertainty in total nitrogen
deposition on sensitive Natura 2000 habitats. Note that the absolute
uncertainties in the uncalibrated OPS model are greater than the
absolute uncertainties on the GDN maps, whereas the relative
uncertainties are approximately equal. This is caused by the fact that
the overestimation by OPS results in a greater denominator in the
computation of the relative uncertainty.
Table 6 Average absolute and relative model uncertainties in total nitrogen
deposition on Natura 2000 habitats.

𝐹𝐹𝑡𝑡𝑡𝑡𝑡𝑡
mol/h/yr
Calibrated
(𝛥𝛥 = 0)
OPS
GDN
4.4

1562
1663
1565

2
�𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚
+ 𝛥𝛥2
𝐹𝐹𝑡𝑡𝑡𝑡𝑡𝑡
mol/h/yr
%
Calibrated;
excluding 𝑠𝑠𝑜𝑜𝑜𝑜𝑜𝑜

2
�𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚
+ 𝛥𝛥2

491

31%

519
491

31%
31%

2
�𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐
+ 𝛥𝛥2
𝐹𝐹𝑡𝑡𝑡𝑡𝑡𝑡
mol/h/yr
%
Calibrated;
including 𝑠𝑠𝑜𝑜𝑜𝑜𝑜𝑜

2
�𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐
+ 𝛥𝛥2

507

32%

534
508

32%
32%

Uncertainty of the national average (total) nitrogen deposition
So far, we have described how the random uncertainties in (total)
nitrogen deposition fluxes on a local scale are obtained. In some
applications (e.g. https://www.clo.nl/en/indicators/en0189-nitrogendeposition), GDN maps are used to compute the national average (total)
nitrogen deposition. In the computation of national averages, the
uncertainties resulting from the calibration will dominate. These
uncertainties can be derived from approximating the calibration process
by a simple subtraction of the mean difference between model and
measurement from the national average model result:

Page 47 of 93

RIVM report 2022-0085

〈𝑋𝑋𝑐𝑐𝑐𝑐𝑐𝑐 〉 = 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 − 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑋𝑋𝑜𝑜𝑜𝑜𝑜𝑜 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁
Eq. 22
where 〈𝑋𝑋𝑜𝑜𝑜𝑜𝑜𝑜 − 𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 is a correction computed by averaging the
difference between model and measurement across all 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜
measurement sites, and 〈 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 averages the calibrated result across all
𝑁𝑁𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔 GDN grid cells. Note that for this simplification of calibration, the
correction is constant for all 𝑁𝑁 grid cells.
We can rewrite Eq. 22 as:
〈𝑋𝑋𝑐𝑐𝑐𝑐𝑐𝑐 〉 = 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 + 〈𝑋𝑋𝑜𝑜𝑜𝑜𝑜𝑜 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 − 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁

Eq. 23

In general, when taking the mean of 𝑁𝑁 draws of random variable 𝑋𝑋, the
variance in that mean is (see Appendix D):
2
𝑠𝑠〈𝑋𝑋〉
= �𝜌𝜌 +

1 − 𝜌𝜌
� 〈𝑠𝑠𝑋𝑋2 〉〈𝑋𝑋〉
𝑁𝑁

Eq. 24

where 𝜌𝜌 is the average correlation between the draws, 〈𝑠𝑠𝑋𝑋2 〉 is the
2
is the uncertainty in the mean of 𝑋𝑋.
uncertainty in 𝑋𝑋, and 𝑠𝑠〈𝑋𝑋〉
Assuming that 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 , 〈𝑋𝑋𝑜𝑜𝑜𝑜𝑜𝑜 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 and 〈𝑋𝑋𝑚𝑚𝑚𝑚𝑚𝑚 〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 are independent, we
can write the uncertainty in 〈𝑋𝑋𝑐𝑐𝑐𝑐𝑐𝑐 〉 as:
2
𝑠𝑠〈𝑋𝑋
𝑐𝑐𝑐𝑐𝑐𝑐 〉

1 − 𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
1 − 𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
2
2
〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 + �𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 +
〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁
� 〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
� 〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
𝑁𝑁𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔
𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜
1 − 𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
2
〉𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁
+ �𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 +
� 〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜
= �𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 +

Eq. 25

The number of GDN grid cells is very large, while the average correlation
coefficient between uncertainties in these grid cells within the
Netherlands is probably low (𝑁𝑁𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔 ≫ 1 − 𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 ), which means that the first
term is negligible. For the observations, it is safe to assume that their
uncertainties are independent (𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 = 0). For the model uncertainties,
we assume that the observation sites are sufficiently far apart to ensure
independence (𝜌𝜌𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋 = 0). The uncertainty in the national average now
reads:
2
𝑠𝑠〈𝑋𝑋
=
𝑐𝑐𝑐𝑐𝑐𝑐 〉

2
2
2
〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
〉 + 〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
〉 〈𝑠𝑠𝑋𝑋𝑋𝑋𝑋𝑋𝑋𝑋
〉
=
𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜
𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜

Eq. 26

where 𝑠𝑠𝑋𝑋2𝑐𝑐𝑐𝑐𝑐𝑐 is the uncertainty in the model plus the uncertainty in the
measurement, as derived for all measured quantities in Chapter 3. For
𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 we take the number of measurement sites against which the
calibration is performed.
Over the 2005-2021 period, the average number of MAN sites that were
used to calibrate the modelled NH3 concentrations against amounted to
184. However, the MAN network itself is calibrated against an average
of 4.6 sites in the LML network. Therefore, we use 𝑁𝑁obs = 4.6 for
NH3 concentrations. The average number of sites that is used to
calibrate the wet deposition fluxes against is 𝑁𝑁obs = 7.9. The
concentrations of NH4 and NO3 are calibrated using an average of
4.1 sites, and the NO2 concentrations are calibrated against 38.3 sites
Page 48 of 93

RIVM report 2022-0085

(excluding the passive samplers). By applying Eq. 26 to the
uncertainties derived in Chapter 3, we obtain the uncertainties in the
averaged calibrated concentrations and wet deposition fluxes listed in
Table 7.
In order to translate the uncertainty in the national average
concentrations to uncertainties in the national average dry deposition
fluxes, we need to account for the fact that on larger spatial scales,
uncertainties in concentrations and dry deposition velocities can no
longer be assumed to be independent, due to mass conservation. For
example, an overestimation in 𝑉𝑉d , will results in lower atmospheric
concentrations and consequently lower dry and wet deposition fluxes.
In Chapter 5, we show that a relative change in 𝑉𝑉d results in a reduced
relative change in the dry deposition fluxes of NHx and NOy. This
reduced sensitivity is expressed by sensitivity factor 𝑓𝑓. A change of 𝑥𝑥%
in the dry deposition velocity yields a change in the dry deposition of 𝑓𝑓 ⋅
𝑥𝑥%. For NHx, the sensitivity factor is approximately 0.38-0.45 and for
NOy, it is approximately 0.63-0.75.
To account for the reduced sensitivity, we multiply the uncertainties in
𝑉𝑉d by these sensitivity factors. Although changes in the dry deposition
velocities also affect the wet deposition fluxes and concentrations, we do
not incorporate these dependencies in our uncertainty analysis. This is
because both wet deposition fluxes and concentrations are calibrated
with measurements and are thus unaffected by biases in 𝑉𝑉d .
In addition to the reduced sensitivity of the modelled deposition to
changes in the dry deposition velocity, we need to account for the fact
that dry deposition fluxes are currently only partially calibrated. That is,
we only correct for biases in the concentrations, while a possible bias in
the dry deposition velocity is not removed from the model result. To
account for the resulting possible bias in the national average dry
deposition flux, we add the uncertainty in the emission totals as derived
by Wever et al. (2023). For NH3, the relative uncertainty in the emission
total is 14% (1 sigma), and for NOx the relative uncertainty is 9%.
Multiplying these relative uncertainties by both domestic and foreign
contributions to the dry deposition fluxes yields uncertainties of 83
mol/ha/yr for NH3 and 32 mol/ha/yr for NOx.
Table 7 presents the resulting uncertainties of the national averages.
The uncertainty in the national average total nitrogen deposition is
201 mol/ha/yr, or 13%. Using the upper and lower bounds of the
sensitivity factors of Vd yields slightly lower and higher uncertainties,
respectively. Hence, we conclude that the uncertainty in the national
average is probably between 10% and 15%.

Page 49 of 93

RIVM report 2022-0085

Table 7 Composition of uncertainties in the national average total nitrogen
deposition. Dry and wet deposition fluxes are denoted as Fd and Fw, and
emission fluxes as FQ. Sensitivity factor f is set to 0.415 for NHx and to 0.69 for
NOy, which are the mid points of the factors derived in Chapter 5.

𝐹𝐹𝑑𝑑,𝑁𝑁𝑁𝑁3
𝐹𝐹𝑑𝑑,𝑁𝑁𝑁𝑁𝑥𝑥
𝐹𝐹𝑑𝑑,𝑁𝑁𝑁𝑁4+
𝐹𝐹𝑤𝑤,𝑁𝑁𝑁𝑁𝑥𝑥
𝐹𝐹𝑤𝑤,𝑁𝑁𝑁𝑁𝑦𝑦

RSD
F
%
26%
39%
30%

SD F
mol/ha/
yr
147
88
9

11%

39%

40%

49

𝐹𝐹𝑐𝑐𝑐𝑐𝑐𝑐
mol/ha/yr
575
227
29

𝐹𝐹𝑂𝑂𝑂𝑂𝑂𝑂
mol/ha/yr
657
239
36

120

106

-15

415
196

406
219

-8
+24

6%
5%
14%
9%

24
9
83
32

1562

1663

147(+3)*

13%

201

𝐹𝐹𝑑𝑑,𝐻𝐻𝐻𝐻𝐻𝐻 +𝑁𝑁𝑁𝑁−
3

𝑅𝑅𝑅𝑅𝑅𝑅 𝐶𝐶
%
11%
2%
18%

f*(RSD
𝑉𝑉𝑑𝑑 )
%
23%
39%
23%

estimated
bias
mol/ha/yr
+82
+11
+7

3

𝐹𝐹𝑄𝑄,𝑁𝑁𝑁𝑁3
𝐹𝐹𝑄𝑄,𝑁𝑁𝑁𝑁𝑥𝑥

sum

*For GDN the average known bias in the national average total nitrogen deposition over
the period 2005-2021 is +3 mol/ha/yr.

4.5

Conclusions
In this chapter, we computed the uncertainty in the total nitrogen
deposition by summing the uncertainties in the individual fluxes derived
in Chapter 3, while accounting for correlations between the uncertainties
of these fluxes. The uncertainty in the total nitrogen deposition depends
on the relative contributions made by the six fluxes at a given location,
but the average uncertainty in the total deposition on sensitive
Natura 2000 habitats is approximately 900-1100 mol/ha/yr (2𝜎𝜎), or
approximately 60-70%. This uncertainty is dominated by uncertainties
in the dry deposition velocities. As a result, removing biases through
calibration or subtracting measurement uncertainties has minor effects
on the uncertainty estimate.
In addition, we have derived an uncertainty in the national average total
nitrogen deposition of 20-30% (2𝜎𝜎) for the calibrated GDN maps. This
uncertainty is dominated by uncertainties resulting from calibration, and
by uncertainties due unknown biases in the dry deposition fluxes. The
uncertainty in the national average total nitrogen deposition is lower
than the local uncertainty within Natura 2000 habitats. This is caused by
negative atmospheric feedbacks, which we will describe in the next
chapter, and by a reduction of the uncertainty due to averaging.

Page 50 of 93

RIVM report 2022-0085

5

Impact of atmospheric feedback on sensitivity of nitrogen
deposition to deposition velocity

Authors: Cor Jacobs, Ronald Hoogerbrugge
5.1

Introduction

5.1.1

Background
Uncertainties in model parameters may only have limited effects on
model estimates for deposition on larger scales. This is because the
atmospheric mass balance implies a negative feedback cycle. For
example, consider a simulation in which the modelled deposition
velocity, 𝑉𝑉d , is too low. In that case, the initial (local) estimate of the
N deposition will also be too low. Too little nitrogen will be removed
from the atmosphere, causing too high N concentrations on large spatial
scales. However, the higher concentrations enhance the deposition
again, compensating the initial effect.
Van Jaarsveld (2004) studied the impact of such surface-atmosphere
interactions on the deposition of NHx, using OPS-Pro 4.1. He reduced 𝑉𝑉d
of ‘grassland’ (pasture) by—roughly – 20%. This land use class
represents almost 40% of the total surface area of the Netherlands.
Thus, a significant influence of surface-atmosphere feedback was
anticipated. Van Jaarsveld (2004) evaluated deposition at six locations:
three locations in an area dominated by ‘grassland’ and three in an area
dominated by cropland. His simulations confirmed that the reduction of
𝑉𝑉d over ‘grassland’ resulted in larger atmospheric concentrations. In
turn, this led to an increase of dry deposition: the effect of the 𝑉𝑉d
change was smaller than anticipated without any change in
concentration. The increased atmospheric concentration also enhanced
the wet deposition. At the ‘grassland’ sites, these effects of the surfaceatmosphere feedback compensated the initial reduction of the dry
deposition by the change in 𝑉𝑉d to a large degree.

Interestingly, the reduced 𝑉𝑉d of ‘grassland’ also resulted in a larger and
opposite effect at the cropland locations. This also attributed to the
afore-mentioned feedback mechanism. At the cropland locations, hardly
any change in 𝑉𝑉d occurred, whereas increased atmospheric
concentrations resulted in enhanced total deposition. This finding by Van
Jaarsveld (2004) reveals another possible implication of surface
atmosphere interactions for uncertainties in N deposition estimates on
larger scales.

In regions with under- or overestimated deposition, negative surfaceatmosphere feedback can reduce the impact on total deposition
estimates, but it cannot completely undo it. In addition, due to
atmospheric transport, N concentrations – and therefore deposition –
may be too high or too low in areas with initially correct. This effect
must be opposite to the effect in the ‘source area’ of the error. As a
consequence, averaging over large areas that include both effects
implies compensation between negative and positive impacts, and hence
a reduction in the total error.
Page 51 of 93

RIVM report 2022-0085

5.1.2

Goal
Here, we assess the effect of the mass balance on error propagation
across scales. The study by Van Jaarsveld (2004) is extended by
examining the sensitivity of the deposition of NHx (NH3 + NH4+) and NOy
(NO + NO2 + NO32- + HNO3-) to 𝑉𝑉d on scales where land-atmosphere
feedback may be significant. We consider effects on atmospheric
concentration and related impacts on both dry and wet deposition and
their sum (total deposition) at 362 receptor points in the Netherlands.
The receptor points represent locations in the MAN and LML
observational networks.

5.2

Method

5.2.1

Model simulations: general
Calculations were performed with OPS-LT version L-5010 for the year
2019. As primary compounds, both NH3 and NO2 were considered.
Emissions for the year 2019 used here include both domestic and
foreign sources of these compounds. Spatially averaged meteorological
conditions in 2019 for the Netherlands were used to drive the model.

5.2.2

Variation of Vd
Deposition velocity 𝑉𝑉d for ‘grassland’ was varied by increasing or
decreasing surface resistance, 𝑅𝑅c . 𝑅𝑅c and 𝑉𝑉d are related by:
𝑉𝑉d =

1
𝑅𝑅a + 𝑅𝑅b + 𝑅𝑅c

Eq. 27

where 𝑅𝑅a is the aerodynamic resistance, determined by turbulent
transport in the lower atmosphere, and 𝑅𝑅b is the quasi-laminar boundary
layer resistance, determined by the depth and conditions of the thin
layer of air neighbouring the surface elements within a canopy (leaves,
soil, water).
OPS uses the DEPAC module (DEPosition of Acidifying Compounds, see
Sauter et al., 2020; Zanten et al., 2010) to compute 𝑅𝑅c . Here, 𝑅𝑅c results
from a weighted contribution by several component resistances: the
canopy-integrated stomatal resistance and the external leaf resistance,
the soil resistance and the in-canopy resistance. Van Jaarsveld (2004)
changed the in-canopy resistance, but in practice, this approach only
results in cases with reduced 𝑉𝑉d . Here, we chose to allow both negative
and positive 𝑉𝑉d changes by forcing 𝑅𝑅c from DEPAC to lower or higher
values.
In order to vary 𝑅𝑅c , and hence 𝑉𝑉d , the value of 𝑅𝑅c computed by DEPAC at
each timestep was systematically multiplied by a change factor, 𝑓𝑓, to
force a difference in 𝑉𝑉d with the reference run in which 𝑓𝑓 = 1 by
definition. For the sensitivity runs, 𝑅𝑅c was varied between 20% and
500% of the value computed by DEPAC (𝑓𝑓 = 0.2, 0.5, 0.8, 1.25, 2.0 and
5.0). Thus, including the reference run, seven runs were performed for
both NH3 and NO2.

Page 52 of 93

RIVM report 2022-0085

5.2.3

Land use-dependent variation of Vd at receptor points
A given relative (fractional) change in 𝑅𝑅c will result in different relative
changes of 𝑉𝑉d . This depends on both the chemical species under
consideration and the land use characteristics near the receptor point
that is being evaluated. The dependence on land use is obtained
because OPS uses a so-called tiling approach to compute 𝑉𝑉d and the
fluxes at receptor points. In this approach, 𝑉𝑉d is first computed for the
nine land use classes implemented in DEPAC (see Table 8). Next, a
weighted average is computed, using the land use fraction of each class
at the receptor point as the weighing factor. As a result, even upon a
fixed change of 𝑉𝑉d of ‘grassland’ only, the tiling approach results in
different effects on the average 𝑉𝑉d at each location. This depends on the
composition of the land use around the receptor point.
Since we apply a 𝑉𝑉d change to ‘grassland’, the resulting 𝑉𝑉d change
‘observed’ at a specific receptor point will mainly depend on the
‘grassland’ fraction surrounding that location. However, Eq. 27 shows
that the relative (fractional) response of 𝑉𝑉d to a change in 𝑅𝑅c also
depends on the values of 𝑅𝑅a and 𝑅𝑅b . Since these resistances also depend
on the surface characteristics, the sensitivity is further modulated by the
land use surrounding the receptor points. Thus, the tiling approach
allows us to effectively sample many different landscape types in our
numerical experiment.
Table 8 Fractional cover of the nine DEPAC land use classes in the Netherlands
(Sauter et al., 2020), in the receptor footprints and the number of footprints
where the specific class is dominant

Land use class

Fractional
cover (%) in
the
Netherlands
Grassland
37.1
Arable land
24.3
Permanent crops (orchards)
1.2
Coniferous forest
4.2
Deciduous forest
6.2
Water
13.2
Built-up area
11.1
Heather and other nature
2.4
Bare soil
0.3

Fractional
cover (%) in
receptor
footprints
31.7
7.5
0.7
3.7
7.6
8.4
4.3
34.4
1.9

Number of
receptor
where
dominant
117
30
4
10
17
21
11
147
5

Table 8 lists the mean characteristics of the landscapes surrounding the
receptor locations in terms of land use classes distinguished by DEPAC.
The land use throughout the Netherlands is presented for comparison.
The land use characteristics near the receptors are determined from the
land use map applied in the OPS simulations. The map has a resolution
of 250 × 250 m2. In turn, this map is based on LGN7 (Hazeu et al.,
2014). The nine DEPAC land use classes in the actual OPS map are
derived from 39 different land use types distinguished in LGN7,
determined at a resolution of 25 × 25 m2 (Sauter et al., 2020). We
assume that land cover in the OPS map square containing a receptor
point roughly represents the main characteristics within the flux

Page 53 of 93

RIVM report 2022-0085

footprint 3 of that receptor. Hence, we will refer to this map area as
‘footprint’.
Most of the receptor points are part of the MAN network. The dominant
land use at a large fraction of the locations is ‘Heather and other
nature’. Some overrepresentation of this land use class in the receptor
sample is justified, because it corresponds to the areas with habitats
sensitive to excessive nitrogen deposition. Cropland appears to be
particularly underrepresented. Especially when average effects on the
national level need to be evaluated, a correction for representativeness
needs to be made.
5.2.4

Analysis of results
At the receptor points, we analyse the annual average of 𝑉𝑉d , the
atmospheric concentration at a height of 3.5 m, and dry and wet
deposition fluxes as well as their sums (total deposition). The present
version of OPS provides the sum of the deposition of primary and
secondary compounds in the output. Therefore, 𝑉𝑉d is evaluated using the
dry deposition of NHx and NOy and the concentration of primary plus
secondary compounds. Thus, because the deposition and concentrations
are lumped, 𝑉𝑉d should be understood as an effective deposition velocity
of NHx and NOy, respectively. In the case of NHx, 𝑉𝑉d also includes a
contribution by the NH3 compensation concentration at the surface.
As we stated before, in the present setup, the impact of changes in 𝑅𝑅c
on 𝑉𝑉d and dry deposition will primarily depend on the fraction of
‘grassland’ area surrounding the receptor point. However, the
composition of the land use beyond the footprint affects the results
through the atmospheric feedback, via its impact on both dry and wet
deposition.

5.3

Results and discussion

5.3.1

Deposition velocity Vd
First, we examine the effect of a change in 𝑅𝑅c (and hence in 𝑉𝑉d ) of
‘grassland’ on the average 𝑉𝑉d perceived at the receptor points, which
have different fractions of ‘grassland’ in their footprint (Figure 10). This
is the starting point of the analysis. The different colours of the symbols
correspond to the different relative changes in 𝑅𝑅c used to drive the 𝑉𝑉d
change over ‘grassland’ (20%,…,500%); the reference case is
represented by the black line (100%). As expected, reduced 𝑅𝑅c results in
increased 𝑉𝑉d and vice versa.
The change in 𝑉𝑉d of homogeneous ‘grassland’, due to the 𝑅𝑅c change
imposed on the numerical experiment, can be estimated from the
results for receptor points with 100% ‘grassland’ in their footprint. For
NHx, the average change in 𝑉𝑉d varies between approximately -20% and
+12% for an increase or decrease of 𝑅𝑅c , respectively, with a factor of 5.
In the case of NOy, the corresponding changes are much larger,
approximately between -35% and +100%. Based on (16), the smaller
response for NHx is expected because of its smaller 𝑅𝑅c relative to 𝑅𝑅a and
𝑅𝑅b . When 𝑅𝑅c is small, 𝑉𝑉d is dominated by 𝑅𝑅a + 𝑅𝑅b . This implies a smaller
(relative) sensitivity to 𝑅𝑅c .
3
When measuring concentrations or fluxes, such as deposition, at a given position, the result contains
contributions from a wider surface area surrounding this position, depending on the measurement height and
meteorological conditions. This area is called the ‘footprint’.

Page 54 of 93

RIVM report 2022-0085

The 𝑉𝑉d ratio at the receptor points varies roughly linearly with the
‘grassland’ fraction in the footprint of the receptor. As explained before,
this is to be expected because of the tiling approach in OPS, based on
weighted contributions to 𝑉𝑉d per land cover type. Deviations from this
linear behaviour are caused by differences in the average 𝑅𝑅a + 𝑅𝑅b , also
related to land use in the footprints of the individual receptors, by
differing concentrations of secondary compounds and, in the case of
NHx, by differences in the NH3 compensation point. The latter quantity
depends on the atmospheric concentration and hence on the large-scale
spatial distribution of sources and sinks. Hence, even for a footprint of
100% ‘grassland’, some variation in 𝑉𝑉d change is obtained.

Figure 10 Change in 𝑉𝑉𝑑𝑑 at the receptor points (expressed as the 𝑉𝑉𝑑𝑑 ratio in the
reference case) as a function of the ‘grassland’ fraction in the footprint of the
receptors. Upper panel: NHx; lower panel: NOy. The purpose of the regression
lines through the results for 𝑅𝑅𝑐𝑐 at 20% and 500% is only to facilitate orientation
and comparison.
Page 55 of 93

Ra�o C at receptor points for NHx

RIVM report 2022-0085

1.15

NHx

1.10
1.05
1.00
0.95
0.90
Rc 20%
Rc 500%

Ra�o C at receptor points for NOy

0.85

1.15 0
1.10

Rc 50%
Rc 200%

Rc 80%
Rc 125%

Reference

10

20

30

40

50

60

70

80

90

100

10

20

30

40

50

60

70

80

90

100

NOyfootprint of the receptor
Grassland frac�on (%) in the
(250m x 250m)

1.05
1.00
0.95
0.90
0.85

0

Grassland frac�on (%) in the footprint of the receptor
(250m x 250m)
Figure 11 Change in atmospheric concentration, C, at the receptor points at a
height of 3.5m (expressed as the ratio with the concentration in the reference
case) as a function of the ‘grassland’ fraction in the footprint of the receptors.
Upper panel: NHx; lower panel: NOy.

5.3.2

Atmospheric concentration
Mass balance considerations imply that the large-scale change in 𝑉𝑉d
must result in corresponding atmospheric concentration changes of
opposite sign (see Section 5.1). This effect is shown in Figure 11, which
depicts the relative change in the atmospheric concentration, C, of
primary plus secondary N compounds at the receptor points, at a height
of 3.5 m. Owing to atmospheric transport and dispersion, the changes
are not restricted to ‘grassland’ areas. For NHx, the concentration
changes typically range between -3% and +5%, with a hardly noticeable
dependence on ‘grassland’ fraction. For NOy, the changes are somewhat
larger on average, corresponding to the larger changes in 𝑉𝑉d , and they
do depend on ‘grassland’ fraction. NOy concentration changes range
between -14% and +7% at 100% ‘grassland’ locations and between 8% and +2% at locations without ‘grassland’.
Page 56 of 93

RIVM report 2022-0085

Figure 12 Change in dry deposition, Ddry, at the receptor points as a function of
the ‘grassland’ fraction in the footprint of the receptors. Ddry is expressed as the
ratio with the dry deposition in the reference case. Upper panel: NHx; lower
panel: NOy

5.3.3

Deposition
The next step in the feedback loop are atmospheric concentration
differences affecting responses of N deposition, both dry (Figure 12) and
wet (Figure 13). At first sight, the relative changes in dry deposition at
the individual receptor points appear to be similar to the relative
changes in 𝑉𝑉d (Figure 10). Indeed, local dry deposition is mainly driven
by local 𝑉𝑉d . However, on closer inspection, the dry deposition ratios turn
out to be smaller than the 𝑉𝑉d ratios. This is caused by the surfaceatmosphere feedback on the larger scales, resulting in the higher
concentrations.

Page 57 of 93

RIVM report 2022-0085

Figure 13 Change in wet deposition, Dwet, at the receptor points as a function of
the ‘grassland’ fraction in the footprint of the receptors. Dwet is expressed as the
relationship with the wet deposition in the reference case. Upper panel: NHx;
lower panel: NOy.

The changes in wet deposition are almost independent of ‘grassland’
fraction, both for NHx and NOy. This is consistent with the idea that
these changes are mainly related large-scale, average concentration
changes.
Finally, the net impact on the total deposition will depend on the
contribution by the component deposition fluxes. Atmospheric feedback
on both dry and wet deposition reduces the relative sensitivity,
depending on the share of dry and wet deposition in the total deposition.

Page 58 of 93

RIVM report 2022-0085

Figure 14 Change in total deposition, Ddry+wet, at the receptor points as a
function of the ‘grassland’ fraction in the footprint of the receptors. Ddry+wet is
expressed as the relationship with the total deposition in the reference case.
Upper panel: NHx; lower panel: NOy.

Figure 14 depicts the change in total deposition on the set of receptor
points considered here. The figure shows that the relative change in
total deposition is much smaller than the change in deposition velocity
𝑉𝑉d and dry deposition (Figure 14 versus Figure 10 for 𝑉𝑉d and versus
Figure 12 for 𝐹𝐹d ).

The results presented in Figure 14 confirm that total N deposition
estimates are comparatively robust, both for NHx and for NOy. This is
because the sensitivity of N deposition to (large-scale) 𝑉𝑉d is limited and
strongly reduced by atmospheric feedback via the concentrations. Note
that the outcome regarding concentrations is often constrained by
observations.
Page 59 of 93

RIVM report 2022-0085

5.3.4

Carry-over effects
Our results confirm the conclusion by Van Jaarsveld (2004) that a given
fractional change in 𝑉𝑉d of ‘grassland’ results in a much smaller fractional
change in dry deposition and in total deposition over ‘grassland’. Our
study also supports the conclusion by Van Jaarsveld (2004) that a largescale change in 𝑉𝑉d of ‘grassland’ may result in opposite changes in total
deposition over other land use classes, without a 𝑉𝑉d change. In the case
of NHx, we obtain positive changes in total deposition for a negative
change in 𝑉𝑉d over ‘grassland’ at receptor points where ‘grassland’
fractions are smaller than ~40%. In the case of NOy, the same applies
to for receptor points with less than ~15% ‘grassland’ in their footprint.
However, according to our results, N deposition responses over non’grassland’ areas are not necessarily stronger than responses over
(nearly) 100% ‘grassland’ areas. On average, the relative decrease of
total deposition at (nearly) 100% ‘grassland’ sites is almost equal to the
relative increase over non-’grassland’ sites. Furthermore, in the case of
NOy the carry-over effect appears to be much smaller for ‘non-grassland’
sites, instead of larger or nearly equal.

5.3.5

Aggregation to larger spatial scales
To conclude, we further explore the consequences of the surfaceatmosphere feedbacks on large scales in the context of spatial
averaging. In this context, the carry-over effect described above may
help to further increase the robustness of N deposition estimates on
ever larger scales.
To this end, we directly examine the sensitivity of dry, wet and total
deposition to 𝑉𝑉d . This is done in three steps (Figure 15):
1. No averaging
The relative change in dry, wet and total deposition is plotted
versus the relative change in 𝑉𝑉d at each individual receptor point
(circles). We do this for the two most extreme changes in 𝑅𝑅c (set
at 20% and 500%, respectively, of its original value) and hence
𝑉𝑉d over ‘grassland’ (cf. Figure 10). Recall that the local change in
𝑉𝑉d depends on the characteristics of the flux footprint of a specific
receptor point. Also recall that 𝑉𝑉d should be understood as an
effective 𝑉𝑉d . Thus, effects of primary and secondary components
are combined, along with effects of the compensation point in the
case of NHy. This explains why the sign of the 𝑉𝑉d change ratios at
receptor points with low ‘grassland’ fractions in their footprint
may be opposite to the sign of the 𝑉𝑉d change of ‘grassland’.
Regarding total deposition (lower panels) at individual receptor
points, we find the relative change to be – only—about one third
of the relative 𝑉𝑉d change for the larger 𝑉𝑉d changes. For example,
a 30% 𝑉𝑉d change at the receptor points induces a change of only
about 10% in total average deposition. For intermediate changes
in 𝑉𝑉d of 5-15%, there may be hardly any effect on the total
deposition. In addition, there are clusters of data that represent
hardly any local change in 𝑉𝑉d . They represent receptors with
hardly any ‘grassland’ in their footprint. At these points, there
may still be a significant change in the total deposition, in a
direction opposite to the change in ‘grassland’ 𝑉𝑉d . Next, it is
explained how changes in both dry and wet deposition contribute
to this negative feedback behaviour.
Page 60 of 93

RIVM report 2022-0085

Relative changes in dry deposition (upper panels) are, roughly,
10-25% lower than the relative changes in 𝑉𝑉d at the receptor
points, depending on the species considered. While they are
obviously strongly related to the ones in 𝑉𝑉d , they are clearly
lower, because of mass balance considerations, resulting in
compensating changes in the atmospheric concentration.
However, especially in the case of NHx (left panels) there is a
cluster of locations with a near-zero change in local 𝑉𝑉d , but, at
the same time, with a significant change in the dry deposition,
consistent with the carry-over effect discussed earlier.
The relative changes in wet deposition (middle panels) are nearly
constant across the entire range of 𝑉𝑉d changes shown here. This
is to be expected, since wet deposition is driven on larger scales,
not by local 𝑉𝑉d at the receptor points. The sign of the changes is
opposite to the one from the change in 𝑉𝑉d , consistent with the
mass balance considerations related to surface-atmosphere
feedback.
2. Averaging over receptor points with a given dominant DEPAC
land use class (diamonds, see also Table 8).
In all cases, the scatter of the individual points can, to some
extent, be considered as a kind of random ‘noise’. This is mainly
due to the differences in land use characteristics in the flux
footprint at a specific receptor point. Hence, when averaging over
a larger number of receptor points, the related uncertainty is
reduced. This is illustrated in the present step. The averages tend
to cluster more towards the ratio 1 for dry as well as for total
deposition. Note that the average for locations with ‘grassland’ as
dominant land use in their footprint is somewhat outside the
cluster of averages for other dominant land use classes. This is to
be expected since the 𝑉𝑉d ‘forcing’ has been applied to ‘grassland’
in the first place.
3. Estimated response of national mean, created by averaging
results over all receptor points, weighed with fractional
occurrence in the Netherlands of the dominant land use in their
footprint (pluses, see also Table 6).
Deviations from the ratio 1 are reduced further. For NHx in
particular (left), negative deviations are compensated by positive
ones at near-zero 𝑉𝑉d change. As a result, the overall mean total
deposition hardly changes (i.e. by less than 0.1%). The overall
relative mean change in dry deposition is only about half the one
in 𝑉𝑉d . For example, in this case, a 6-7% decrease in mean 𝑉𝑉d
results in a 2-3% decrease in dry deposition. In the case of NOy
(right), the compensation of negative ratios by positive ones is
lower, but still significant. For example, a 28% increase in mean
𝑉𝑉d changes the overall mean dry deposition by about 18% and
the total deposition by about 8%.

Page 61 of 93

RIVM report 2022-0085

Figure 15 Effect of a change in 𝑉𝑉𝑑𝑑 (x-axis) on dry deposition (upper), wet
deposition (middle) and total deposition (lower) for NHx (left) and NOy (right).
Recall that changes are given as weighed, relative changes of Vd and the
deposition considering land use composition in the footprint of a receptor point.
The relative changes are expressed as the ratio of the value after a change in Rc
of pasture ‘grassland’ with a factor of five to the value obtained in the reference
case without a change in ‘grassland’ Rc. Open circles (o): at individual receptor
points; diamonds (): averaged over clusters of receptors with a specific
dominant land use; pluses (+): overall mean. Orange: 𝑅𝑅𝑐𝑐 for ‘grassland’ is 20%
of base value; Blue: 𝑅𝑅𝑐𝑐 for ‘grassland’ is 500% of base value. Note the
differences in axis scales between left (NHx) and right (NOy). The ratio of the
relative change in the average Fd over the relative change in the average Vd
(indicated by the + marks in the upper panels) was used to obtain the sensitivity
factors applied in Section 4.4 to derive the uncertainty in the national average
total nitrogen deposition.
Page 62 of 93

RIVM report 2022-0085

5.4

Conclusions
Total N deposition calculations by dispersion models such as OPS appear
to be comparatively robust on large spatial scales. In the sensitivity
study presented here, large-scale systematic changes in deposition
velocity over ‘grassland’ (~40% of the Dutch land area) are found to to
have hardly any impact on total N deposition estimates on the national
scale. This has two main reasons.
First, mass balance considerations dictate negative feedback between
the surface and the atmosphere: more removal by deposition implies
lower atmospheric concentrations, and vice versa. These changes
modulate both wet and dry deposition, such that they compensate the
initial change. As a result, in our sensitivity studies, the relative change
in total N deposition at specific receptor points was found to be –
typically—only about one third of the relative 𝑉𝑉d change.
Second, owing to atmospheric transport, atmospheric concentration is
affected beyond ‘grassland’ areas as well. Thus, N deposition over non’grassland’ areas will also change. The sign of the change in such areas
tends to be opposite to the one over ‘grassland’. This reduces the effect
of the initial change on the average on a national scale, since negative
and positive deviations in N deposition tend to compensate each other.
Calibration of the model outcomes against measurements can become
confounded by the atmospheric feedbacks. However, on a national
scale, this confounding effect will be limited since the average
calibration effect is very small.

Page 63 of 93

RIVM report 2022-0085

Page 64 of 93

RIVM report 2022-0085

6

Spatial variability of nitrogen deposition: comparing of
GCN/GDN and DN

Authors: Sebastiaan Hazelhorst, Wouter Marra
6.1

Introduction
In this chapter, we compare nitrogen depositions and ammonia
concentrations at two different spatial resolutions. The goal of this
comparison is to gain insight into the effect of spatial resolution of
emissions and grid size on calculated concentrations and depositions.
For this comparison, we use data from DN 2021 and GCN/GDN 2021 of
ammonia concentrations (unpublished for DN) and total nitrogen
deposition. DN contains nitrogen depositions calculated at a highresolution (1-hectare) hexagonal grid covering Dutch Natura 2000
areas. GCN/GDN contains depositions and concentrations on a 1 x 1 km2
grid covering the whole of the Netherlands.
DN and GCN/GDN use the same model (OPS) and emission data to
calculate nitrogen depositions. The difference between the two sets is
the resolution of the input and the output. For the calculation of GDNGCN, gridded emissions at a resolution of 1 × 1 km2 are used. DN
makes use of more detailed information concerning the locations of
emission sources, where available. Examples of more detailed emission
sources are the locations of roads and the emissions from farms and
agricultural fields.
The higher resolution of the DN maps is used in combination with highresolution data on designated nitrogen-sensitive habitats and breeding
and resting sites for rare and threatened species.

6.2

Methods
For the comparison between DN 2021 and GCN/GDN 2021, we used
several datasets. First, we compared modelled NH3 concentrations with
observations from MAN. For this comparison, we used calculations with
emission files that were also used for DN 2021 and GCN/GDN 2021.
These emission files are elaborated in more detail below. The receptor
file used for the calculations contained the exact coordinates of the
measuring sites, including height. These calculations were performed for
the year 2018, using emissions from the Dutch Emissieregistratie
(Wever et al., 2020), and meteorological conditions for that year. Note
that the modelled concentrations of NH3 used here are not biascorrected (not calibrated against the measurements), as the case with
depositions in DN or GCN/GDN, and therefore directly indicate the
performance of the model.
Furthermore, we also used NH3 concentrations and total nitrogen
depositions from DN 2021 and GCN-GDN 2021 for the projection year
2030. 4 These datasets are obtained using projected emissions. To this
end, data from the Dutch Emissieregistratie is combined with projections
of emissions in the Netherlands from the KEV-2020 (Klimaat en Energie
4
A comparison between AERIUS and GCN/GDN for an actual year was not possible. The reason is that all
calculations for DN are based on long-term average meteorology (2005-2014) and are bias-corrected on the
basis of five years of differences between model and observations.

Page 65 of 93

RIVM report 2022-0085

Verkenning; Climate and Energy Outlook). For other European countries,
emissions from the MACC III emission inventory 5 are combined with
projections of the NAPCP scenario from the EU’s second Clean Air
Outlook (IIASA, 2021) 6.
The calculations for the projection year are performed using average
weather conditions from the 2005-2014 period. The results were biascorrected on the basis of five years of differences between model results
and measurements from the 2014-2018 period. More details on the
construction of the DN and the GCN/GDN deposition maps can be found
in the AERIUS manual (https://www.aerius.nl/nl/handboeken) and the
GCN/GDN report (Hoogerbrugge et al., 2022), respectively. In the
analyses, we only used data on so-called relevant hexagons. Relevant
hexagons are, in short, hexagons containing nitrogen-sensitive habitats
and breeding and resting sites for rare and threatened species.
In order to link DN data with GCN/GDN data, a geographic overlay was
made between the hexagonal grid structure of DN and the square grid of
GCN/GDN. For every GCN/GDN grid cell, the mean total nitrogen
deposition of DN hexagons within each grid cell is calculated.
Thus, we compare the DN values with those of GCN/GDN in the same
grid cell. The analysis only includes GDN grid cells containing more than
90 relevant hexagons (>90% coverage). GCN/GDN grid cells with fewer
than ninety relevant hexagons will generally not give a fair comparison,
as these hexagons are often located on one side of a GCN/GDN grid cell
only. Therefore, such clusters of hexagons cannot be considered
representative for the GCN/GDN grid cell they are located in. This
analysis essentially shows the difference between calculating the
deposition at a coarse resolution and calculating the deposition at a high
resolution and averaging those values for a larger area.
In addition to a comparison of averages, we directly compared total
nitrogen depositions and ammonia concentrations between individual
relevant hexagons and the corresponding GCN/GDN grid cells. This
provides an assessment of the sub-grid spatial variability of depositions
and concentrations within GCN/GDN grid cells. This variability is mainly
related to the effect of local emission sources and small-scale (100 m)
variations in terrain properties on the outcome. The differences between
the two resolutions are further illustrated with an example for the
Natura 2000 area Meijendal & Berkheide, located in the province of
Zuid-Holland.
6.3

Results

6.3.1

Comparing of modelled NH3 concentrations with observations
Modelled concentrations calculated in DN 2021 and GCN/GDN 2021 are
each lower than the observations (Figure 16). With a linear regression
coefficient of 0.91, the GCN/GDN setup performs slightly better than
DN, which has a coefficient of 0.85.
Differences between these two setups can be explained by the locations
of emissions. As GCN/GDN uses emissions aggregated on a scale of
1 × 1 km2 and DN a more detailed emissions on a hectare scale, the
location of the measuring site within a GCN/GDN grid cell plays a key
role here (Figure 17).
5
6

https://cordis.europa.eu/project/id/633080/reporting
https://ec.europa.eu/environment/air/clean_air/outlook.htm

Page 66 of 93

RIVM report 2022-0085

Figure 16 Modelled (y-axis) and measured (x-axis) NH3 concentrations for the
year 2018. In red, concentrations modelled by DN 2021. In blue, concentrations
are from GCN/GDN 2021. Observations are yearly averages from MAN and LML.
Modelled data are not bias-corrected. The trendlines are based on linear
interpolation with the formula 𝑦𝑦 = 𝑎𝑎 ⋅ 𝑥𝑥. The standard deviations of a in the
formula are equal to 0.015 for DN and 0.016 for GCN.

When, for instance, the measuring site is relatively close (within
approximately a couple of hundred meters) to an emission location, the
concentration modelled by DN is higher than the one modelled with
GCN. This is because the same amount of emission is spread over a
large surface in GCN, whereas in DN, the NH3 is emitted at a specific
location. In contrast, as the distance between emission location and
measurement site increases, the modelled concentration from GCN/GDN
becomes relatively higher than in DN, since the modelled concentration
from GCN/GDN remains roughly the same, while in DN, it decreases
significantly with distance.
The dominant wind direction in relation to the location of the local
sources and receptor also plays a role in explaining the differences
between DN and GCN/GDN. At measurement sites that are located on a
large distance from sources, i.e. in a grid cell without large emission
sources, modelled concentrations of DN and GCN/GDN are more
comparable. For the calibrations of both maps, measurement sites that
are located directly next to large sources are excluded.

Page 67 of 93

RIVM report 2022-0085

Figure 17 Example of the effect of more detailed emission locations on the ratio
between GCN/GDN and DN modelled concentrations on measuring sites (in
yellow the location of Binnenveld is shown here). In red, the outlines of the
1 × 1 km2 GCN/GDN grid cells are shown. Emissions for GCN/GDN are spread
equally over these grid cells. The pink circles are the emission sources used in
the DN setup.

6.3.2

Comparing of DN hexagonal grid and GDN square grid.
Deposition data of 2030 from DN compares very well with GDN
(regression coefficient of 1.003 and R2 of 0.9996, Figure 19). A root
mean square error of 25 mol N/ha/yr is relatively small compared to an
average deposition of approximately 1220 mol/ha/yr. The largest
differences are found at locations with transitions in terrain properties
(Figure 19). For example, on the island Terschelling, the effect of large
differences in surface roughness due to the land-sea transition is clearly
visible in the average deposition.
Although the average total nitrogen deposition of DN hexagons in a
1 × 1 km2 grid cell are often within 5% (Figure 19), local variations in
deposition are larger (Figure 20). As a result of local terrain properties
and emission sources, more detailed patterns in deposition become
visible. For example, the local forest area close to the beach area
Wassenaarse Slag has a higher deposition in DN than in GDN. This is a
result of the larger local roughness length in the more detailed
calculations in DN.

Page 68 of 93

RIVM report 2022-0085

Figure 18 Comparison of total N deposition between GDN 2021 and DN 2021 for
the year 2030 for grid cells containing more than ninety relevant hexagons. The
trendline is based on linear interpolation with the formula 𝑦𝑦 = 𝑎𝑎 ⋅ 𝑥𝑥. The standard
deviation of a in the formula is equal to 0.00061.

Figure 19 Relative difference in N deposition between GDN 2021 and DN 2021
for the year 2030, for grid cells containing more than ninety relevant hexagons.
Page 69 of 93

RIVM report 2022-0085

Figure 20 Total nitrogen deposition in 2030 on the square grid of 1 × 1 km2 of
GDN (left) and on the DN hexagonal grid of 1 hectare (middle) for the coastal
Natura 2000 area Meijendel & Berkheide. An aerial picture is shown in the right
picture as a reference. Although the average differences between both grids are
small, more detailed patterns in deposition due to local terrain properties and
emission sources become visible in the DN deposition map.

Figure 21 Total nitrogen deposition in 2030 calculated on DN hexagons of
1 hectare (y-axis) compared to the deposition of the corresponding GDN grid cell
of 1 × 1 km2 (x-axis). The locations in the Natura2000 area ‘Meijendal &
Berkheide’ are highlighted in red. The data is filtered on relevant hexagons and
grid cells containing more than ninety hexagons. The trendline is based on linear
interpolation with the formula 𝑦𝑦 = 𝑎𝑎 ⋅ 𝑥𝑥. The standard deviation of a in the
formula is equal to 0.00043.

When comparing depositions on 1-hectare hexagons directly to the GDN
grid cell it is located in, the local variability within a 1 × 1 km2 grid cell
can be quantified (Figure 21). A root mean square error of
178 mol N/ha/yr corresponds to an average variability in nitrogen
deposition of 14.6% when compared to the average deposition
(~ 1220 mol/ha/yr on the selection of GDN grid cells). Furthermore, the
RMSE of 178 mol/ha/yr indicates that local variation in deposition can be
Page 70 of 93

RIVM report 2022-0085

quite large. This is especially important for small Natura 2000 areas, as
local variations can have a large effect on the average deposition of such
a small area. The variability is attributed to a combination of more
detailed emission locations and local terrain properties.

Figure 22 NH3 concentration in 2030 calculated on DN hexagons of 1 hectare (yaxis) compared to the NH3 concentration of the corresponding GCN/GDN grid cell
of 1 × 1 km2 (x-axis). The locations in the Natura2000 area Meijendal &
Berkheide are highlighted in red. The data is filtered on relevant hexagons and
grid cells containing more than ninety hexagons. The trendline is based on linear
interpolation with the formula 𝑦𝑦 = 𝑎𝑎 ⋅ 𝑥𝑥. The standard deviation of a in the
formula is equal to 0.00025.

For a large part, the variability in total nitrogen deposition in Figure 21
is caused by variations in deposition velocity due to varying terrain
properties. However, part of the variability is a result of differences in
emissions. In order to separate the effect of terrain properties and that
of emissions, we also looked at concentrations on different grid levels
(Figure 22). A root mean square error of 0.26 ug m-3 NH3 indicates an
average variability in NH3 concentration of 8.3% when compared to the
average concentration (~ 3.12 ug m-3 NH3 on the selection of GDN grid
cells). This means that the variation in the deposition in a GDN grid is
largely caused by the variation in the terrain properties: two thirds of
the variation in the deposition in a GCN/GDN grid cell is explained by the
variation in the deposition calculation that considers more detailed
terrain properties.

Page 71 of 93

RIVM report 2022-0085

6.4

Conclusions
In general, GDN and DN have no systematic difference for either the
calculated NH3 concentrations or the nitrogen deposition. Variations in
the 1 × 1 km2 scale of the GDN grid can be significant. An RMSE
between 1 ha hexagons within the 1 × 1 km2 grid average of nearly
10% is observed for NH3 concentrations, and of approximately 15% for
nitrogen deposition. The major part of this variation originates from
terrain properties.

Page 72 of 93

RIVM report 2022-0085

7

General discussion and conclusion

This status report combines several approaches to studying the
uncertainty in the determination of the nitrogen deposition in the
Netherlands. Compared to the uncertainty analyses of calculated
concentrations, the analysis of deposition is much more complicated.
Firstly, because several paths of deposition exist of which six
compounds are considered in this report. Secondly, because the
parameter of interest, nitrogen deposition, is not measured directly for
three of the compounds. For these, the uncertainty in deposition is
determined by combining an uncertainty in the concentration estimate
with an uncertainty in the deposition velocity.
In combining the uncertainties of the various compounds, another
complexity arises because some estimates should be made on the
correlation between the compounds. Since the compounds are
calculated using the same model, OPS 5.1.0.1, the assumption that
uncertainties are uncorrelated is not realistic. For example, the
important equation for the deposition resistance has similar 𝑅𝑅a for all
airborne compounds. A proxy of the correlation in uncertainties in
deposition velocities has been estimated using the correlation between
the spatial distribution of the compounds. If we use this estimate, the
uncertainty estimate increases from 23 to 32%. Obtaining better
knowledge of this correlation is therefore recommended.
The dry depositions of NH3, NO2, NO3-/NH4+, have high relative
uncertainties due to a large uncertainty in the deposition velocity.
The most important component is dry deposition of NH3. For NH3, a
large set of concentration measurements is available, combined with a
smaller set of deposition measurements. The standard deviation of the
differences between the measured and modelled deposition is quite
large. The result indicates a systematic difference between the modelled
and measured depositions. However, we assess the information as being
insufficient for the determination of the bias and subsequent calibration.
This assessment of insufficient coherent information was based on
several observations:
• additional measurements in the Netherlands, for shorter periods,
show opposite differences;
• measurement values from international studies show a wide
range of results, even after selection of similar land use classes
and correction for concentration differences.
A more detailed analysis of the differences between the calculated and
measured dry deposition at the COTAG sites will be carried out.
Without the possibility for a bias correction, a large uncertainty in
NH3 deposition and deposition velocity is estimated, resulting in an
overall uncertainty estimate of 60-70% (2𝜎𝜎) on a local scale.
In the estimate of the uncertainty of the total N deposition, the
contribution of the NH3 dry deposition is dominating. Fortunately, more
deposition measurement locations have been planned. Both the
concentration measurements and the deposition measurements are
concentrated in Natura 2000 areas. This implies that the uncertainty
Page 73 of 93

RIVM report 2022-0085

estimates are representative for locations like the Natura 2000 areas.
Uncertainties for other locations, such as agricultural areas, built up
areas and transition locations, might be different.
In this report, six compounds are considered. In the OPS model, some
compounds actually contain several sub-compounds, which may behave
quite differently. For example, deposition of secondary oxidised nitrogen
is considered to consist of the relatively slow deposition of ammonium
nitrate particles and the quick deposition of HNO3. It is recommended to
adjust OPS to distinguish between these compounds.
The above considerations all reflect the small-scale uncertainty
estimates and their correlations. In Chapter 5, atmospheric feedbacks
are studied that may affect uncertainty on large scales. This is achieved
by introducing variations (‘deviations’) in the deposition velocity over
‘grassland’. In general, the variation in the deposition velocity over
‘grasslands’ results in a negative feedback. This can be explained as
follows: decreasing the deposition velocity will decrease the deposition,
which will result in increased concentration. The higher concentration
then leads to a higher deposition and thus partly compensates the initial
effect. As a result, the decrease of deposition at ‘grassland’ locations is
less than might be expected on the basis of the change in the deposition
velocity.
At locations with no or very small amounts of ‘grasslands’, only an
increase in concentration will be observed, leading to somewhat
increased deposition. As a result, the average deposition over larger
areas (consisting of ‘grassland’ as well as other land cover) will be even
less sensitive to a change in the deposition velocity than expected by
the local effects.
Finally, a compensation through the wet deposition is calculated as well.
The feedback part of the study shows that the impact of the uncertainty
in the deposition velocity on the total national nitrogen deposition is
considerably smaller than the estimate on the local scale. On the basis
of this feedback study, the uncertainty in the total national nitrogen
deposition is estimated at 20-30%.

Page 74 of 93

RIVM report 2022-0085

8

References

Bohrnstedt, G.W. and Goldberger, A.S., 1969. On the exact covariance
of products of random variables. Journal of the American
Statistical Association, 64(328): 1439-1442,
https://doi.org/10.2307/2286081.
Cape, J.N. et al., 2008. Estimate of annual NH3 dry deposition to a
fumigated ombrotrophic bog using concentration-dependent
deposition velocities. Atmospheric Environment, 42(27): 66376646, https://doi.org/10.1016/j.atmosenv.2008.04.027.
Dobben, H.v. and Hinsberg, A.v., 2008. Overzicht van kritische
depositiewaarden voor stikstof, toegepast op habitattypen en
Natura 2000-gebieden Alterra, Wageningen,
https://edepot.wur.nl/45419.
EMEP, 2022. Transboundary particulate matter, photo-oxidants,
acidification and eutrophication components. , Joint MSC-W &
CCC & CEIP &CIAM Report. EMEP Status Report 1/2022,
https://emep.int/publ/reports/2022/EMEP_Status_Report_1_202
2.pdf.
Flechard, C.R. et al., 2011. Dry deposition of reactive nitrogen to
European ecosystems: a comparison of inferential models across
the NitroEurope network. Atmospheric Chemistry and Physics,
11(6): 2703-2728, https://doi.org/10.5194/acp-11-2703-2011.
Goodman, L.A., 1960. On the exact variance of products. Journal of the
American Statistical Association, 55(292): 708-713,
https://doi.org/10.2307/2281592.
Hazeu, G.W., Schuiling, C., Dorland, G.J., Roerink, G.J., Naeff, H.S.D.
and Smidt, R.A., 2014. Landelijk Grondgebruiksbestand
Nederland versie 7 (LGN7); Vervaardiging, nauwkeurigheid en
gebruik, Alterra Wageningen UR (University & Research centre),
Wageningen.
Hoogerbrugge, R. et al., 2022. Grootschalige concentratie- en
depositiekaarten Nederland. Rapportage 2022,
https://doi.org/10.21945/RIVM-2022-0059.
Hordijk, L. et al., 2020. Meer meten, robuuster rekenen - Eindrapport
van het Adviescollege Meten en Berekenen Stikstof,
https://www.rijksoverheid.nl/documenten/rapporten/2020/06/15
/meer-meten-robuuster-rekenen.
Marra, W., Hazelhorst, S., Brandt, K., Kruit, R.W. and Schram, J., 2022.
Monitor stikstofdepositie in Natura 2000-gebieden 2022.
Uitgangssituatie voor de Wet Stikstofreductie en
Natuurverbetering, Rijksinstituut voor Volksgezondheid en Milieu,
Bilthoven, https://doi.org/10.21945/RIVM-2022-0120.
Nguyen, P. and Wesseling, J., 2016. Metingen van
stikstofdioxideconcentraties (NO2) met Palmes buisjes - Periode
2012-2015, Rijksinstituut voor Volksgezondheid en Milieu,
Bilthoven, https://www.rivm.nl/bibliotheek/rapporten/20160089.pdf.

Page 75 of 93

RIVM report 2022-0085

Noordijk, H., Braam, M., Rutledge-Jonker, S., Hoogerbrugge, R., Stolk,
A.P. and van Pul, W.A.J., 2020. Performance of the MAN
ammonia monitoring network in the Netherlands. Atmospheric
Environment, 228: 117400,
https://doi.org/10.1016/j.atmosenv.2020.117400.
Rutledge-Jonker, S. et al., 2023. Ammonia deposition measured with
Conditional Time-Averaged Gradient (COTAG) systems in the
Netherlands National Institute for Public Health and the
Environment (RIVM), Bilthoven, https://doi.org/10.21945/RIVM2022-0202.
Sauter, F., Sterk, M., Swaluw, E.v.d., Kruit, R.W., Vries, W.d. and Pul,
A.v., 2020. The OPS-model; Description of OPS 5.0.0.0, National
Institute for Public Health and the Environment (RIVM),
Bilthoven, https://www.rivm.nl/documenten/uitgebreidemodelbeschrijving-van-ops-versie-5000.
Sauter, F., Sterk, M., Swaluw, E.v.d., Kruit, R.W., Vries, W.d. and Pul,
A.v., 2023. The OPS-model; Description of OPS 5.1.1.0, National
Institute for Public Health and the Environment (RIVM),
Bilthoven, https://www.rivm.nl/documenten/uitgebreidemodelbeschrijving-van-ops-versie-5110.
Schrader, F. and Brümmer, C., 2014. Land use specific ammonia
deposition velocities: a review of recent studies (2004-2013).
Water Air Soil Pollut, 225(10): 2114,
https://doi.org/10.1007/s11270-014-2114-7.
Siteur, K., Volten, H., Vendel, K., Bongers, P. and Hoogerbrugge, R.,
2021. Nauwkeurige NO2-concentratiekaarten dankzij vrijwilligers.
Tijdschrift Lucht 17(3): 12-16.
Sutton, M.A., Dragosits, U., Geels, C., Gyldenkærne, S., Misselbrook,
T.H. and Bussink, W., 2015. Review on the scientific
underpinning of calculation of ammonia emission and deposition
in the Netherlands,
https://www.tweedekamer.nl/kamerstukken/detail?id=2015D324
94&did=2015D32494.
van Jaarsveld, J.A., 2004. The Operational Priority Substances model.
Description and validation of OPS-Pro 4.1, Rijksinstituut voor
Volksgezondheid en Milieu, Bilthoven,
https://www.rivm.nl/bibliotheek/rapporten/500045001.pdf.
Velders, G.J.M., Aben, J.M.M., Jaarsveld, J.A.v., Pul, W.A.J.v., Vries,
W.J.d. and Zanten, M.C.v., 2010. Grootschalige stikstofdepositie
in Nederland. Herkomst en ontwikkeling in de tijd, Planbureau
voor de Leefomgeving (PBL), Den Haag/Bilthoven,
https://www.rivm.nl/bibliotheek/rapporten/500088007.pdf.
Vendel, K.J.A. et al., 2023. Dry deposition of ammonia in a coastal dune
area: Measurements and modeling. Atmospheric Environment,
https://doi.org/10.1016/j.atmosenv.2023.119596: 119596,
https://doi.org/10.1016/j.atmosenv.2023.119596.
Walker, J.T. et al., 2019. Aspects of uncertainty in total reactive
nitrogen deposition estimates for North American critical load
applications. Science of The Total Environment, 690: 1005-1018,
https://doi.org/10.1016/j.scitotenv.2019.06.337.
Wever, D. et al., 2020. Informative Inventory Report 2020 : Emissions
of transboundary air pollutants in the Netherlands 1990-2018,
Rijkinstituut voor Volksgezondheid en Milieu, Bilthoven,
https://doi.org/10.21945/RIVM-2020-0032.
Page 76 of 93

RIVM report 2022-0085

Wever, D. et al., 2023. Informative Inventory Report 2023. Emissions of
transboundary air pollutants in the Netherlands 1990–2021,
Rijkinstituut voor Volksgezondheid en Milieu, Bilthoven,
https://doi.org/10.21945/RIVM-2023-0057.
Wichink Kruit, R. et al., 2021. Op weg naar een optimale meetstrategie
voor stikstof, Rijksinstituut voor Volksgezondheid en Milieu,
Bilthoven, https://doi.org/10.21945/RIVM-2021-0118.
Wichink Kruit, R., Braam, M., Hoogerbrugge, R. and Pul, A.v., 2020.
Implementation of a data fusion approach to assess the
concentration and dry deposition of ammonia in the Netherlands,
Rijksinstituut voor Volksgezondheid en Milieu, Bilthoven,
https://doi.org/10.21945/RIVM-2020-0076.
Wintjen, P., Schrader, F., Schaap, M., Beudert, B. and Brümmer, C.,
2022a. Forest–atmosphere exchange of reactive nitrogen in a
remote region – Part I: Measuring temporal dynamics.
Biogeosciences, 19(2): 389-413, https://doi.org/10.5194/bg-19389-2022.
Zanten, M.C.v., Sauter, F.J., Kruit, R.J.W., Jaarsveld, J.A.v. and Pul,
W.A.J.v., 2010. Description of the DEPAC module. Dry deposition
modelling with DEPAC_GCN2010, Rijksinstituut voor
Volksgezondheid en Milieu, Bilthoven,
https://www.rivm.nl/bibliotheek/rapporten/680180001.pdf.

Page 77 of 93

RIVM report 2022-0085

Page 78 of 93

RIVM report 2022-0085

Appendix A Proportional random uncertainty and alternative
models

In our analysis, we estimate the uncertainty in the concentrations and
deposition fluxes at locations other than the observation sites by
assuming proportional uncertainty:
𝑠𝑠𝑋𝑋2 (𝑋𝑋) = RSD2 𝑋𝑋 2

Eq. A1

where 𝑋𝑋 is a measured quantity, 𝑠𝑠𝑋𝑋 is the model uncertainty of that
quantity, and RSD is the relative uncertainty. We approximate parameter
𝑠𝑠𝑋𝑋
. In the literature (Thunis, 2013; Pernigotti 2013), the
RSD by RSD = 〈𝑥𝑥〉
uncertainty is often split into a part that is constant and a part that is
proportional to 𝑋𝑋:
𝑠𝑠𝑋𝑋2 (𝑋𝑋) = 𝑠𝑠02 + RSD2 𝑋𝑋 2

Eq. A2

In this appendix, we study the validity of Eq. A1, as used in the present
analysis.
By subdividing variable 𝑋𝑋 into twenty equally-sized bins, with a
minimum of five observations per bin, and plotting both the Eq. A1 and
the uncertainties per bin, we obtain Figure A1. From visual inspection,
the proportional model used in our study seems to perform well for
relevant ranges of concentrations and wet deposition fluxes.
Besides a visual comparison, we computed the Akaike Information
Criterion (AIC) for both models Eq. A1 and Eq. A2. The model with the
lowest AIC gives the most likely relation between 𝑠𝑠𝑋𝑋 and 𝑋𝑋. For this
analysis, we estimate the parameters of Eq. A2 through linear
regression.
We find that the proportional model (Eq. A1) gives a lower AIC for all
compounds except NH3. For NH3, the AIC of the more advanced model
of Eq. A2 is lower. For simplicity and reproducibility, we use Eq. A1 in
this report for all compounds.

Page 79 of 93

RIVM report 2022-0085

Figure A1 Binned uncertainties (points) and the proportional model of Eq. A1
(lines). In black: model uncertainty + measurement uncertainty (𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐 ) In red:
only model uncertainty (𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 ).

Page 80 of 93

RIVM report 2022-0085

Appendix B Independence between deposition velocities and
concentrations

In Chapters 3 and 4, we assume independence between the residuals in
deposition velocities and concentrations. To test this assumption,
deposition velocities have been computed using NH3 concentration
measurements at ~5m above the surface at the COTAG sites. The
concentration measurements are part of the MAN network and were
performed independently from measurements by the COTAGs. We found
no correlation between residuals of the calibrated model concentrations
and the residuals of the modelled deposition velocities (𝜌𝜌 = −0.006). This
is illustrated by Figure A2.

Figure A2 Residuals in the effective dry deposition velocity and calibrated
NH3 concentration at the COTAG sites.

Note that although the covariance between residuals in 𝐶𝐶 and 𝑉𝑉𝑑𝑑 may be
negligible on a local scale (i.e. the footprint of the COTAG), this may not
hold true on larger spatial scales. On larger scales, atmospheric
feedbacks between deposition and concentration, involving depletion
due to deposition, become important. This is studied in detail in
Chapter 5.

Page 81 of 93

RIVM report 2022-0085

Appendix C Uncertainties in dry deposition fluxes and
correlations between fluxes

C.1

Variance of the product of two independent random variables
In many cases, we need to estimate the variance of combinations of
variables, where we know the variances of these underlying variables.
For the sum of variables, the resulting variance can be determined in a
simple way. Determining the products of two independent random
variables is more complex, but an exact solution is available. If two
variables 𝑥𝑥 and 𝑦𝑦 are independent, the variance of their product is
(Goodman, 1960):
Var(𝑥𝑥𝑥𝑥) = 𝐸𝐸(𝑥𝑥)2 Var(𝑦𝑦) + 𝐸𝐸(𝑦𝑦)2 Var(𝑥𝑥) + Var(𝑥𝑥)Var(𝑦𝑦)

Eq. A3

𝑠𝑠𝐹𝐹2𝑑𝑑 = 𝐶𝐶 2 𝑠𝑠𝑉𝑉2𝑑𝑑 + 𝑉𝑉𝑑𝑑2 𝑠𝑠𝐶𝐶2 + 𝑠𝑠𝑉𝑉2𝑑𝑑 𝑠𝑠𝐶𝐶2

Eq. A4

substituting 𝑉𝑉𝑑𝑑 , 𝐶𝐶 and 𝐹𝐹𝑑𝑑 (= 𝑉𝑉𝑑𝑑 ∙ 𝐶𝐶) the uncertainty in 𝐹𝐹𝑑𝑑 can be written
as:

or as in the main text:

C.2

𝑠𝑠𝐹𝐹2𝑑𝑑 = 𝐹𝐹𝑑𝑑2 �

𝑠𝑠𝑉𝑉2𝑑𝑑

2 +

𝑉𝑉𝑑𝑑

𝑠𝑠𝐶𝐶2 𝑠𝑠𝑉𝑉2𝑑𝑑 𝑠𝑠𝐶𝐶2
+
�
𝐶𝐶 2 𝑉𝑉𝑑𝑑2 𝐶𝐶 2

Eq. A5

Exact Covariance of Products of Random Variables
In 1969, Bohrnstedt and Goldberger published an analysis of the
covariance of products of random variables. Assuming multivariate
normality, they showed that:
Cov(𝑥𝑥𝑥𝑥, 𝑢𝑢𝑢𝑢)
= 𝐸𝐸(𝑥𝑥)𝐸𝐸(𝑢𝑢)Cov(𝑦𝑦, 𝑣𝑣) + 𝐸𝐸(𝑥𝑥)𝐸𝐸(𝑣𝑣)Cov(𝑦𝑦, 𝑢𝑢) + 𝐸𝐸(𝑦𝑦)𝐸𝐸(𝑢𝑢)Cov(𝑥𝑥, 𝑣𝑣)
+ 𝐸𝐸(𝑦𝑦)𝐸𝐸(𝑣𝑣)Cov(𝑥𝑥, 𝑢𝑢) + Cov(𝑥𝑥, 𝑢𝑢)Cov(𝑦𝑦, 𝑣𝑣) + Cov(𝑥𝑥, 𝑣𝑣)Cov(𝑦𝑦, 𝑢𝑢)

Eq. A6

Where 𝑥𝑥, 𝑦𝑦, 𝑢𝑢 and 𝑣𝑣 are random variables.
This result can be used to determine the covariance between two dry
deposition fluxes from the underlying deposition velocities and
concentrations. Substituting 𝐶𝐶, 𝑉𝑉𝑑𝑑 and 𝐹𝐹𝑑𝑑 and rewriting yields:
Cov�𝐹𝐹𝑑𝑑,𝑖𝑖 , 𝐹𝐹𝑑𝑑,𝑗𝑗 � = 𝐹𝐹𝑑𝑑,𝑖𝑖 𝐹𝐹𝑑𝑑,𝑗𝑗 �

Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 ,𝑉𝑉𝑑𝑑,𝑗𝑗 �
𝑉𝑉𝑑𝑑,𝑖𝑖 𝑉𝑉𝑑𝑑,𝑗𝑗

+

Cov�𝐶𝐶𝑖𝑖 ,𝐶𝐶𝑗𝑗 �
𝐶𝐶𝑖𝑖 𝐶𝐶𝑗𝑗

+

Cov�𝐶𝐶𝑖𝑖 ,𝑉𝑉𝑑𝑑,𝑗𝑗 �
𝐶𝐶𝑖𝑖 𝑉𝑉𝑑𝑑,𝑗𝑗

+

Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 ,𝐶𝐶𝑗𝑗 �
𝑉𝑉𝑑𝑑,𝑖𝑖 𝐶𝐶𝑗𝑗

Cov�𝐶𝐶𝑖𝑖 ,𝐶𝐶𝑗𝑗 �Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 ,𝑉𝑉𝑑𝑑,𝑗𝑗 �+Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 ,𝐶𝐶𝑗𝑗 �Cov�𝐶𝐶𝑖𝑖 ,𝑉𝑉𝑑𝑑,𝑗𝑗 �
𝐹𝐹𝑑𝑑,𝑖𝑖 𝐹𝐹𝑑𝑑,𝑗𝑗

�

+

Eq. A7

Assuming 𝐶𝐶 and 𝑉𝑉𝑑𝑑 can be treated as independent, this simplifies to:
Cov�𝐹𝐹𝑑𝑑,𝑖𝑖 , 𝐹𝐹𝑑𝑑,𝑗𝑗 �
Cov�𝐶𝐶𝑖𝑖 , 𝐶𝐶𝑗𝑗 � Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 , 𝑉𝑉𝑑𝑑,𝑗𝑗 � Cov�𝐶𝐶𝑖𝑖 , 𝐶𝐶𝑗𝑗 �Cov�𝑉𝑉𝑑𝑑,𝑖𝑖 , 𝑉𝑉𝑑𝑑,𝑗𝑗 �
= 𝐹𝐹𝑑𝑑,𝑖𝑖 𝐹𝐹𝑑𝑑,𝑗𝑗 �
+
+
�
𝐶𝐶𝑖𝑖 𝐶𝐶𝑗𝑗
𝑉𝑉𝑑𝑑,𝑖𝑖 𝑉𝑉𝑑𝑑,𝑗𝑗
𝐹𝐹𝑑𝑑,𝑖𝑖 𝐹𝐹𝑑𝑑,𝑗𝑗

Page 82 of 93

Eq. A8

RIVM report 2022-0085

or in terms of correlation:
𝜌𝜌�𝐹𝐹𝑑𝑑,𝑖𝑖 , 𝐹𝐹𝑑𝑑,𝑗𝑗 � =

𝐹𝐹𝑑𝑑,𝑖𝑖 𝐹𝐹𝑑𝑑,𝑗𝑗
𝑠𝑠𝐹𝐹𝑑𝑑,𝑖𝑖 𝑠𝑠𝐹𝐹𝑑𝑑,𝑖𝑖

𝑠𝑠𝑉𝑉𝑑𝑑,𝑖𝑖 𝑠𝑠𝑉𝑉𝑑𝑑,𝑗𝑗
𝑠𝑠𝐶𝐶𝑖𝑖 𝑠𝑠𝐶𝐶𝑗𝑗
+ 𝜌𝜌�𝑉𝑉𝑑𝑑,𝑖𝑖 , 𝑉𝑉𝑑𝑑,𝑗𝑗 �
𝐶𝐶𝑖𝑖 𝐶𝐶𝑗𝑗
𝑉𝑉𝑑𝑑,𝑖𝑖 𝑉𝑉𝑑𝑑,𝑗𝑗
𝑠𝑠𝐶𝐶𝑖𝑖 𝑠𝑠𝐶𝐶𝑗𝑗
𝑠𝑠𝑉𝑉𝑑𝑑,𝑖𝑖 𝑠𝑠𝑉𝑉𝑑𝑑,𝑗𝑗
+ 𝜌𝜌�𝐶𝐶𝑖𝑖 , 𝐶𝐶𝑗𝑗 �
𝜌𝜌�𝑉𝑉𝑑𝑑,𝑖𝑖 , 𝑉𝑉𝑑𝑑,𝑗𝑗 �
�
𝐶𝐶𝑖𝑖 𝐶𝐶𝑗𝑗
𝑉𝑉𝑑𝑑,𝑖𝑖 𝑉𝑉𝑑𝑑,𝑗𝑗

�𝜌𝜌�𝐶𝐶𝑖𝑖 , 𝐶𝐶𝑗𝑗 �

Eq. A9

Note that the correlation coefficient between two dry deposition fluxes is
solely determined by the relative uncertainties in 𝐶𝐶, 𝑉𝑉𝑑𝑑 and 𝐹𝐹𝑑𝑑 , as
derived in Chapter 3, and by the correlations between the
concentrations and the dry deposition velocities. These correlation
coefficients are determined in Appendix C.3.
From the result of Bohrnstedt and Goldberger (1969), the correlation
between a wet and a dry deposition flux can also be derived:
𝜌𝜌�𝐹𝐹𝑤𝑤,𝑖𝑖 , 𝐹𝐹𝑑𝑑,𝑗𝑗 � =

𝑠𝑠𝐶𝐶𝑗𝑗
𝑠𝑠𝑉𝑉𝑑𝑑,𝑗𝑗
𝐹𝐹𝑤𝑤,𝑖𝑖
�𝜌𝜌�𝐹𝐹𝑤𝑤,𝑖𝑖 , 𝐶𝐶𝑗𝑗 �
+ 𝜌𝜌�𝐹𝐹𝑤𝑤,𝑖𝑖 , 𝑉𝑉𝑑𝑑,𝑗𝑗 �
�
𝑠𝑠𝐹𝐹𝑤𝑤,𝑖𝑖
𝐶𝐶𝑗𝑗
𝑉𝑉𝑑𝑑,𝑗𝑗

𝜌𝜌�𝐹𝐹𝑤𝑤,𝑖𝑖 , 𝐹𝐹𝑑𝑑,𝑗𝑗 � =

𝑠𝑠𝐶𝐶𝑗𝑗
𝐹𝐹𝑤𝑤,𝑗𝑗
𝜌𝜌�𝐹𝐹𝑤𝑤,𝑖𝑖 , 𝐶𝐶𝑗𝑗 �
𝑠𝑠𝐹𝐹𝑤𝑤,𝑖𝑖
𝐶𝐶𝑗𝑗

Eq. A10

Assuming dry deposition velocities and dry deposition fluxes to be
independent yields:

C.3

Eq. A11

Correlation between residuals in concentrations, dry deposition
velocities and wet deposition fluxes
Correlations between residuals (modelled minus observed value) can be
derived for pairs of measured quantities, using sites where both are
measured over the same period. Table A1 shows the number of such
overlapping observations for the 2005-2021 period.
Note that the number of overlapping observations is generally much
smaller than the number of observations for each individual compound
(diagonal in Table A1).
The correlations between residuals of measured variables are presented
in Table A2.

Page 83 of 93

RIVM report 2022-0085

Table A1 Number of observation pairs that can be used to compute the
correlation coefficient between residuals of the measured quantities.

𝑪𝑪𝐍𝐍𝐍𝐍𝟑𝟑

𝒏𝒏

𝑪𝑪𝐍𝐍𝐍𝐍𝟑𝟑

3123

𝑪𝑪𝐍𝐍𝐍𝐍𝟐𝟐

𝑪𝑪𝐍𝐍𝐍𝐍𝟐𝟐

155

𝑪𝑪𝐍𝐍𝐍𝐍𝟒𝟒+

43

𝑪𝑪𝐍𝐍𝐍𝐍−𝟑𝟑

43

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱

39

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

793

67

67

71

71

70

70

52

52

70

52

52

135

135

𝑪𝑪𝐍𝐍𝐍𝐍𝟒𝟒+

𝑪𝑪𝐍𝐍𝐍𝐍−𝟑𝟑

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

39

135

Table A2 Correlation coefficients of the residuals, for each combination of the
measured variables

𝝆𝝆

𝑪𝑪𝐍𝐍𝐍𝐍𝟑𝟑

𝑪𝑪𝐍𝐍𝐍𝐍𝟑𝟑

1.00

𝑪𝑪𝐍𝐍𝐍𝐍𝟐𝟐

𝑪𝑪𝐍𝐍𝐍𝐍𝟐𝟐

0.17

𝑪𝑪𝐍𝐍𝐍𝐍𝟒𝟒+

0.38

𝑪𝑪𝐍𝐍𝐍𝐍−𝟑𝟑

-0.12

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱

-0.08

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

1.00

0.03

-0.17

-0.07

-0.30

1.00

0.77

0.56

0.40

1.00

0.69

0.51

1.00

0.67

𝑪𝑪𝐍𝐍𝐍𝐍𝟒𝟒+
𝑪𝑪𝐍𝐍𝐍𝐍−𝟑𝟑

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

-0.27

1.00

For the approach outlined in Appendix C.2, correlation coefficients
between the residuals of dry deposition velocities are required. As the
dry deposition velocities are currently not measured, they cannot be
computed. There is much overlap in the model variables that determine
the dry deposition velocity, and therefore the errors in the model are
probably highly correlated. Here, we approximate the required
correlation coefficients by computing the correlation of modelled
deposition velocities on 1 x 1 km2 GDN grid cells over the 2005-2021
period, weighted by the fraction of relevant Natura 2000 habitats. This
results in Table A3.
Table A3 Correlation coefficients of the residuals, for each combination of the dry
deposition velocities

𝝆𝝆

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑
1.00

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟐𝟐

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+
𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍−𝟑𝟑

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟐𝟐
0.81

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+
0.71

𝑽𝑽𝒅𝒅,𝐍𝐍𝐍𝐍−𝟑𝟑

1.00

0.74

0.67

1.00

0.79

0.85

1.00

Combining both tables, using the equations derived in Appendix C.2,
yields two correlation tables. One that is based on the combined model
and observation uncertainty (𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐 , Table A4) and one that is solely based
Page 84 of 93

RIVM report 2022-0085

on the model uncertainty (e.g. after subtracting the observation
uncertainty, 𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 , Table A5). Note that the correlation coefficients in
both tables are roughly equal.
Table A4 Correlation coefficients of the residuals, for each combination of
deposition fluxes. The computation of this table is described in Appendix C.2,
using relative uncertainties in 𝐹𝐹𝑑𝑑 , 𝐹𝐹𝑤𝑤 , 𝑉𝑉𝑑𝑑 and C, before removing the
measurement uncertainty (𝑠𝑠𝑐𝑐𝑐𝑐𝑐𝑐 ).

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑

𝝆𝝆

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

1.00

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱

0.73
1.00

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−
𝟑𝟑

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

0.58
0.54
1.00

𝟑𝟑

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

1.00

0.30
1.00

0.22
0.67
1.00

𝟑𝟑

0.66
0.55
0.76

𝟑𝟑

-0.03
-0.02
0.34

-0.10
-0.07
0.24

Table A5 Correlation coefficients of the residuals, for each combination of the
deposition fluxes. The computation of this table is described in Appendix C.2,
using relative uncertainties in 𝐹𝐹𝑑𝑑 , 𝐹𝐹𝑤𝑤 , 𝑉𝑉𝑑𝑑 and C, after removing the measurement
uncertainty (𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚 ).

𝝆𝝆

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟑𝟑

1.00

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−
𝟑𝟑

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱
𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

𝟑𝟑

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝐱𝐱

0.74
1.00

𝑭𝑭𝒅𝒅,𝐍𝐍𝐍𝐍𝟒𝟒+

0.58
0.55
1.00

𝑭𝑭𝒅𝒅,𝐇𝐇𝐇𝐇𝐇𝐇 +𝐍𝐍𝐍𝐍−

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐱𝐱

𝑭𝑭𝒘𝒘,𝐍𝐍𝐍𝐍𝐲𝐲

1.00

0.30
1.00

0.22
0.67
1.00

𝟑𝟑

0.66
0.57
0.76

𝟑𝟑

-0.03
-0.01
0.34

-0.10
-0.05
0.24

Page 85 of 93

RIVM report 2022-0085

Appendix D The average of multiple draws of a random
variable

The average of 𝑁𝑁 draws of a random variable is:
〈𝑋𝑋〉𝑁𝑁 =

𝑁𝑁

1
� 𝑋𝑋𝑖𝑖
𝑁𝑁

Eq. A12

𝑖𝑖=1

We are interested in the variance of the average of 𝑁𝑁 random variables,
which can be written as:
𝑁𝑁

𝑁𝑁

𝑁𝑁

𝑖𝑖=1

𝑖𝑖=1 𝑗𝑗=1

1
1
Var(〈𝑋𝑋〉𝑁𝑁 ) = 2 Var �� 𝑋𝑋𝑖𝑖 � = 2 � � Cov(𝑋𝑋𝑖𝑖 , 𝑋𝑋𝑗𝑗 )
𝑁𝑁
𝑁𝑁

Eq. A13

Where Cov() is the covariance between two draws of 𝑋𝑋. The number of
terms on the diagonal of the covariance matrix is 𝑁𝑁. The total
contribution by these terms to the variance in the average is 𝑁𝑁 ∙ Var(𝑋𝑋).
The number of cross terms is 𝑁𝑁 2 − 𝑁𝑁. The total contribution by the cross
terms is (𝑁𝑁 2 − 𝑁𝑁)𝜌𝜌Var(𝑋𝑋), where 𝜌𝜌 is the average correlation between
draws of 𝑋𝑋. Hence, the variance in the average of 𝑁𝑁 random variables is:
Var(〈𝑋𝑋〉𝑁𝑁 ) =

𝑁𝑁

𝑁𝑁

1
� � Cov�𝑋𝑋𝑖𝑖 , 𝑋𝑋𝑗𝑗 �
𝑁𝑁 2
𝑖𝑖=1 𝑗𝑗=1

1
1 − 𝜌𝜌
= 2 �𝑁𝑁 ∙ Var(𝑋𝑋) + (𝑁𝑁 2 − 𝑁𝑁)𝜌𝜌Var(𝑋𝑋)� = �𝜌𝜌 +
� Var(𝑋𝑋)
𝑁𝑁
𝑁𝑁

Page 86 of 93

Eq. A14

RIVM report 2022-0085

Appendix E

E.1

Assessing Vd of NH3 in an international context
Figure A3 shows a comparison between 𝑉𝑉𝑑𝑑 values obtained from the
international scientific literature for vegetation surfaces ((semi)natural,
agricultural, forest) (Schrader and Brümmer, 2014; henceforth SB14),
from Dutch measurement sites and from OPS runs for the Dutch sites.
From the SB14 review, only data derived from estimates in all four
seasons per year is shown. Most of the 𝑉𝑉𝑑𝑑 data in their set was obtained
with the so-called inferential method. This is considered a modelling
approach, and therefore, data from this method is distinguished from
data based on micrometeorological measurements (open versus closed
symbols). Little information on 𝑉𝑉𝑑𝑑 from more direct observations is
available, particularly for ‘forest’.
It is clear that within the land use classes considered here, the spread in
𝑉𝑉𝑑𝑑 values listed by SB14 is considerable. Note that some of the
variability will be caused by the fact that differences in measurement
height and climatic conditions cannot be accounted for. Nonetheless, in
the ‘agriculture and seminatural’ class, the variability for the subset of
measured data (filled symbols) remains similar to the one including the
inferential method. For ‘forest’, only three comparatively high 𝑉𝑉𝑑𝑑 values
remain when only considering measurements.
Observations from the Dutch sites show less spread. This may be due to
the fact that climatic conditions and heights for the evaluation of 𝑉𝑉𝑑𝑑 are
similar, with the exception of the Solleveld site. Regarding OPS results,
we conclude that the range in 𝑉𝑉𝑑𝑑 from the simulations agrees reasonably
well within the range of both Dutch and international observations.
Again, the Solleveld results may be an exception, since OPS has
overestimated 𝑉𝑉𝑑𝑑 for that site. Possible causes are discussed in Vendel et
al. (2023). So far, no direct comparison can be made between OPS and
measurements at forest sites. Details on the data from SB14 and from
Vendel et al. (2023) can be found in Appendix E.2.

Page 87 of 93

RIVM report 2022-0085

Figure A3 Comparison of measured 𝑉𝑉𝑑𝑑 values from the international review by
Schrader and Brümmer (2014; SB14) with Dutch measurements at COTAG sites
(Rutledge-Jonker et al., 2023) and the Solleveld site (Vendel et al., 2023)
(labelled together as NL OBS) and with corresponding OPS computations for the
Dutch sites (OPS). These data sources have been highlighted in the top of the
figure, above the corresponding 𝑉𝑉𝑑𝑑 data clusters. Diamonds: SB14 seminatural;
squares: SB14 agricultural; circles: SB14 forest (from left to right: coniferous,
deciduous, mixed). Open symbols refer to 𝑉𝑉𝑑𝑑 values derived from the inferential
method, which is essentially a modelling approach; coloured symbols refer to 𝑉𝑉𝑑𝑑
values derived from micrometeorological measurements. X: Dutch
measurements (orange: Bargerveen and Oostelijke Vechtplassen (n=10); green:
Solleveld (n=1); purple: Hoge Veluwe (n=2)); +: output from OPS with colour
coding as for measurements. See Appendix E.2 for further details.

E.2

Interpretation and processing of data in Schrader & Brümmer
(2014) and Vendel et al. (2023)
To further assess the uncertainty of deposition velocity, 𝑉𝑉𝑑𝑑 values for
NH3 collected by SB14 were used. In their review, these authors
extracted 𝑉𝑉𝑑𝑑 data published in the scientific literature in the 2004-2013
period. In total, they collected 61 𝑉𝑉𝑑𝑑 values from 42. In 59 cases, the 𝑉𝑉𝑑𝑑
value could be assigned to one of the seven land use classes they
distinguished. The land use classes are presented in Table A6, along
with the number of 𝑉𝑉𝑑𝑑 values reported in each class. Individual 𝑉𝑉𝑑𝑑 data
and metadata from the review has been made available as
Page 88 of 93

RIVM report 2022-0085

Supplementary Information to the paper, partly to encourage additional
analyses.
Table A6 Land use classes distinguished by Schrader and Brümmer (2014) and
number of observations in each class upon filtering.

Land use

Agricultural
Coniferous forest
Deciduous forest
Mixed forest
Seminatural
Urban
Water

All
(n=59)

𝒘𝒘𝒔𝒔 = 𝟒𝟒
(n=44)

18
6
4
4
19
5
3

10
5
3
4
15
5
2

𝒘𝒘𝒔𝒔 = 𝟒𝟒 &
measured
(n=16)
1
1
1
1
7
3
2

The 𝑉𝑉𝑑𝑑 data listed in SB14 is meant to be used as annual averages.
However, in practice, some 𝑉𝑉𝑑𝑑 values are based on a limited number of
observations throughout the year. To this end, SB14 also provide quality
flags indicative of the data coverage within a year, with two weight
indicators (𝑤𝑤):
• The 𝑤𝑤𝑠𝑠 indicator weighs the number of seasons that were
sampled: (1-4).
Subscript 𝑠𝑠 refers to ‘season’.
• The 𝑤𝑤𝑐𝑐 indicator refers to the measurement duration within the
sampled season.
From 1 for ‘1 day to 3 weeks per season’ to 4 for ‘nine weeks or
more per season’.
Subscript 𝑐𝑐 refers to ‘data coverage’.

Since our purpose is a comparison of annual 𝑉𝑉𝑑𝑑 , we decided to exclude
data with 𝑤𝑤𝑠𝑠 < 4. This choice implies that 𝑉𝑉𝑑𝑑 should have been obtained
from observations in all seasons of the year. It can be argued that, as
long as each season has been sampled, the exact duration of the
measurement periods in each of the seasons is less important. Thus, 𝑤𝑤𝑐𝑐
was not used to further filter the data.

In our own analysis, we also decided to consider the method on which
the 𝑉𝑉𝑑𝑑 was based. First, data from earlier reviews was excluded because
it may contain unwanted and unknown overlaps among the data
considered. Second, much of the remaining 𝑉𝑉𝑑𝑑 data in SB14 is based on
the so-called inferential method (INF). Measurements—e.g.
concentrations, standard meteorological conditions and sometimes also
turbulence observations—are used to drive a 𝑉𝑉𝑑𝑑 model. In fact, this
makes the INF method a modelling approach. Nevertheless, reports on
𝑉𝑉𝑑𝑑 derived from INF have been accepted in the peer-reviewed scientific
literature as defensible assessments of 𝑉𝑉𝑑𝑑 (e.g. Flechard et al., 2011).
Therefore, we decided not to exclude this data beforehand. Instead, we
distinguish it from 𝑉𝑉𝑑𝑑 measurements that are based solely on
micrometeorological measurements, such as eddy covariance
measurements and gradient techniques.

Page 89 of 93

RIVM report 2022-0085

A limited literature search was performed to assess whether new annual
𝑉𝑉𝑑𝑑 measurements have been published since SB14. This was performed
by means of a Google Scholar search of peer-reviewed scientific
publications in English that cite SB14. The idea is that new 𝑉𝑉𝑑𝑑 values will
probably be compared with the numbers collected by SB14.
The search resulted in thirty publications (Table A7; last update: 12
April 2023). Most of these papers cite SB14 in a more general sense, for
example to underpin or illustrate the statement that 𝑉𝑉𝑑𝑑 is uncertain.
Other papers choose a 𝑉𝑉𝑑𝑑 value partly or entirely based on SB14 in order
to use this value to estimate the deposition from their own measured or
computed concentrations. Papers citing SB14 and publishing new 𝑉𝑉𝑑𝑑
values on the basis of measurements are rare. Amongst these, one
paper (Wintjen et al., 2022a) would meet the criterion 𝑤𝑤𝑠𝑠 = 4, but this
paper reported lumped 𝑉𝑉𝑑𝑑 values for total nitrogen, so that 𝑉𝑉𝑑𝑑 of NH3
ammonia could not be separated out anymore. All in all, no data was
added to the SB14 data.
Several shortcomings in the SB14 review can be identified. Major ones
are the lack of climatic information (wind speed), more specific
information on terrain properties (e.g. aerodynamic roughness) and
information on the height for which 𝑉𝑉𝑑𝑑 was calculated. For the time
being, we assume that this contributes to the scatter in the 𝑉𝑉𝑑𝑑 numbers
from SB14.
In summary, it was decided to use from SB14:
1. all land-use specific 𝑉𝑉𝑑𝑑 data, for which 𝑤𝑤𝑠𝑠 = 4;
2. A subset of 1), in which 𝑉𝑉𝑑𝑑 was obtained from
micrometeorological observations.
Figure A3 in Appendix E.1 compares data from SB14 for predominantly
vegetated land use (i.e., excluding the land use classes ‘urban’ and
‘water’) with observations from the COTAG sites (Rutledge-Jonker et al.,
2023) and with measurements of 𝑉𝑉𝑑𝑑 near Solleveld (Vendel et al., 2023).
Results from OPS simulations for the Dutch measurement sites are
compared with the observations for the following years: Bargerveen
(COTAG) 2015-2019; Oostelijke Vechtplassen (COTAG) 2015-2019;
Hoge Veluwe (COTAG) 2018-2019; Solleveld (GRAHAM) 2014/2015.
Note that in the latter case, data from 2014 and 2015 was combined
into one full year. In OPS, all measurement sites are classified as ‘other
nature’ and as such they are compared with the ‘seminatural’ from
SB14. So far, no recent Dutch NH3 observations for any of the forest
types distinguished in SB14 are available. The receptor height in OPS
was set to 5m, which matches the upper measurement level of the
COTAG sites. However, at Solleveld, 𝑉𝑉𝑑𝑑 has been evaluated at the height
of ~1m (Vendel et al., 2023). Thus, it should be kept in mind that 𝑉𝑉𝑑𝑑
from OPS would become somewhat larger when placing the receptor
point on this lower level (since the aerodynamic resistance becomes
smaller).

Page 90 of 93

RIVM report 2022-0085

Papers citing Schrader and Brümmer (2014) published in the English
peer-reviewed scientific literature [obtained with Google Scholar; last
update 12-4-2023]
Adon, M., Galy‐Lacaux, C., Serça, D., Guedant, P., Vongkhamsao, A.,
Rode, W., . . . Guérin, F. (2018). First assessment of inorganic
nitrogen deposition budget following the impoundment of a
subtropical hydroelectric reservoir (Nam Theun 2, Lao PDR).
Journal of Geophysical Research: Atmospheres, 123, 12,413412,428.
Ahrends, B., Schmitz, A., Prescher, A.-K., Wehberg, J., Geupel, M.,
Andreae, H., & Meesenburg, H. (2020). Comparison of methods for
the estimation of total inorganic nitrogen deposition to forests in
Germany. Frontiers in Forests and Global Change, 3, 103.
Aksoyoglu, S., & Prévôt, A. S. (2018). Modelling nitrogen deposition: dry
deposition velocities on various land-use types in Switzerland.
International Journal of Environment and Pollution, 64, 230-245.
Baker, J., Battye, W. H., Robarge, W., Arya, S. P., & Aneja, V. P.
(2020). Modeling and measurements of ammonia from poultry
operations: Their emissions, transport, and deposition in the
Chesapeake Bay. Science of The Total Environment, 706, 135290.
Battye, W. H., Bray, C. D., Aneja, V. P., Tong, D., Lee, P., & Tang, Y.
(2019). Evaluating ammonia (NH3) predictions in the NOAA NAQFC
for Eastern North Carolina using ground level and satellite
measurements. Journal of Geophysical Research: Atmospheres,
124, 8242-8259.
Braun, S., Ahrends, B., Alonso, R., Augustin, S., García-Gómez, H.,
Hůnová, I., . . . Thimonier, A. (2022). Nitrogen deposition in
forests: Statistical modeling of total deposition from throughfall
loads. Frontiers in Forests and Global Change, 5.
Geupel, M., Heldstab, J., Schäppi, B., Reutimann, J., Bach, M.,
Häußermann, U., . . . Breuer, L. (2021). A national nitrogen target
for Germany. Sustainability, 13, 1121.
Guo, H., Liu, J., Froyd, K. D., Roberts, J. M., Veres, P. R., Hayes, P. L., .
. . Weber, R. J. (2017). Fine particle pH and gas–particle phase
partitioning of inorganic species in Pasadena, California, during the
2010 CalNex campaign. Atmospheric Chemistry and Physics, 17,
5703-5719.
Guo, H., Otjes, R., Schlag, P., Kiendler-Scharr, A., Nenes, A., & Weber,
R. J. (2018). Effectiveness of ammonia reduction on control of fine
particle nitrate. Atmospheric Chemistry and Physics, 18, 1224112256.
Guo, H., Sullivan, A. P., Campuzano‐Jost, P., Schroder, J. C., Lopez‐
Hilfiker, F. D., Dibb, J. E., . . . Nenes, A. (2016). Fine particle pH
and the partitioning of nitric acid during winter in the northeastern
United States. Journal of Geophysical Research: Atmospheres,
121, 10,355-310,376.
Guo, X., Pan, D., Daly, R. W., Chen, X., Walker, J. T., Tao, L., . . .
Zondlo, M. A. (2022). Spatial heterogeneity of ammonia fluxes in a
deciduous forest and adjacent grassland. Agricultural and Forest
Meteorology, 326, 109128.

Page 91 of 93

RIVM report 2022-0085

Häni, C., Flechard, C., Neftel, A., Sintermann, J., & Kupper, T. (2018).
Accounting for field-scale dry deposition in backward Lagrangian
stochastic dispersion modelling of NH3 emissions. Atmosphere, 9,
146.
Hůnová, I., Kurfürst, P., Vlček, O., Stráník, V., Stoklasova, P.,
Schovankova, J., & Srbová, D. (2016). Towards a better spatial
quantification of nitrogen deposition: A case study for Czech
forests. Environmental Pollution, 213, 1028-1041.
Jia, Y., Yu, G., Gao, Y., He, N., Wang, Q., Jiao, C., & Zuo, Y. (2016).
Global inorganic nitrogen dry deposition inferred from ground-and
space-based measurements. Scientific Reports, 6, 1-11.
Johnson, J., Cummins, T., & Aherne, J. (2016). Critical loads and
nitrogen availability under deposition and harvest scenarios for
conifer forests in Ireland. Science of The Total Environment, 541,
319-328.
Khan, M. A. H., Lowe, D., Derwent, R., Foulds, A., Chhantyal-Pun, R.,
McFiggans, G., . . . Shallcross, D. (2020). Global and regional
model simulations of atmospheric ammonia. Atmospheric
Research, 234, 104702.
Kosonen, Z., Schnyder, E., Hiltbrunner, E., Thimonier, A., Schmitt, M.,
Seitler, E., & Thöni, L. (2019). Current atmospheric nitrogen
deposition still exceeds critical loads for sensitive, semi-natural
ecosystems in Switzerland. Atmospheric Environment, 211, 214225.
Kubota, T., Kuroda, H., Watanabe, M., Takahashi, A., Nakazato, R.,
Tarui, M., . . . Ouchi, T. (2020). Role of advection in atmospheric
ammonia: A case study at a Japanese lake basin influenced by
agricultural ammonia sources. Atmospheric Environment, 243,
117856.
Lassman, W., Collett Jr, J. L., Ham, J. M., Yalin, A. P., Shonkwiler, K. B.,
& Pierce, J. R. (2020). Exploring new methods of estimating
deposition using atmospheric concentration measurements: a
modeling case study of ammonia downwind of a feedlot.
Agricultural and Forest Meteorology, 290, 107989.
Thimonier, A., Kosonen, Z., Braun, S., Rihm, B., Schleppi, P., Schmitt,
M., . . . Thöni, L. (2019). Total deposition of nitrogen in Swiss
forests: Comparison of assessment methods and evaluation of
changes over two decades. Atmospheric Environment, 198, 335350.
Tiwari, A., Kumar, P., Baldauf, R., Zhang, K. M., Pilla, F., Di Sabatino,
S., . . . Pulvirenti, B. (2019). Considerations for evaluating green
infrastructure impacts in microscale and macroscale air pollution
dispersion models. Science of The Total Environment, 672, 410426.
Van Der Graaf, S., Dammers, E., Segers, A., Kranenburg, R., Schaap,
M., Shephard, M. W., & Erisman, J. W. (2022). Data assimilation of
CrIS NH 3 satellite observations for improving spatiotemporal NH 3
distributions in LOTOS-EUROS. Atmospheric Chemistry and
Physics, 22, 951-972.
Van Der Graaf, S. C., Dammers, E., Schaap, M., & Erisman, J. W.
(2018). How are NH 3 dry deposition estimates affected by
combining the LOTOS-EUROS model with IASI-NH 3 satellite
observations? Atmospheric Chemistry and Physics, 18, 1317313196.
Page 92 of 93

RIVM report 2022-0085

Van Der Graaf, S. C., Kranenburg, R., Segers, A. J., Schaap, M., &
Erisman, J. W. (2020). Satellite-derived leaf area index and
roughness length information for surface–atmosphere exchange
modelling: a case study for reactive nitrogen deposition in northwestern Europe using LOTOS-EUROS v2. 0. Geoscientific Model
Development, 13, 2451-2474.
Walker, J. T., Beachley, G., Zhang, L., Benedict, K. B., Sive, B. C., &
Schwede, D. B. (2020). A review of measurements of air-surface
exchange of reactive nitrogen in natural ecosystems across North
America. Science of The Total Environment, 698, 133975.
Weber, R. J., Guo, H., Russell, A. G., & Nenes, A. (2016). High aerosol
acidity despite declining atmospheric sulfate concentrations over
the past 15 years. Nature Geoscience, 9, 282-285.
Wintjen, P., Schrader, F., Schaap, M., Beudert, B., & Brümmer, C.
(2022a). Forest–atmosphere exchange of reactive nitrogen in a
remote region–Part I: Measuring temporal dynamics.
Biogeosciences, 19, 389-413.
Wintjen, P., Schrader, F., Schaap, M., Beudert, B., Kranenburg, R., &
Brümmer, C. (2022b). Forest–atmosphere exchange of reactive
nitrogen in a remote region–Part II: Modeling annual budgets.
Biogeosciences, 19, 5287-5311.
Yi, W., Shen, J., Liu, G., Wang, J., Yu, L., Li, Y., . . . Wu, J. (2021). High
NH3 deposition in the environs of a commercial fattening pig farm
in central south China. Environmental Research Letters, 16,
125007.
Zöll, U., Brümmer, C., Schrader, F., Ammann, C., Ibrom, A., Flechard,
C. R., . . . Kutsch, W. L. (2016). Surface–atmosphere exchange of
ammonia over peatland using QCL-based eddy-covariance
measurements and inferential modeling. Atmospheric Chemistry
and Physics, 16, 11283-11299.

Page 93 of 93

Published by:
National Institute for Public Health
and the Enviroment, RIVM
P.O. Box 1 | 3720 BA Bilthoven
The Netherlands
www.rivm.nl/en
December 2023

Committed to
health and sustainability