ECN-C--04-027
METHODE ONTWIKKELING VOOR VERBETERING VAN DE VERIFICATIE VAN NIETCO2 BROEIKASGASEMISSIES Een NOVEM-ROB2002 project A.T. Vermeulen A. Hensen A.C.B den Ouden G. Pieterse
DECEMBER 2004
Revisies A B Opgesteld door:
Goedgekeurd door:
A.T. Vermeulen Geverifieerd door:
G.J. de Groot Vrijgegeven door:
A. Bleeker
J.W. Erisman
ECN-Brandstoffen Conversie & Milieu Werkeenheid:
Verantwoording Dit rapport is het eindresultaat van de ECN bijdrage aan het gelijknamige project uitgevoerd in opdracht van Senter/Novem in het kader van het ROB (Reductie Overige Broeikasgassen) programma, projectnummer 0377-02-01-01-008. Het ECN projectnummer was 7.2817.
Abstract The main purpose of the project was the improvement of the accuracy of emission verification strategies for non-CO2 greenhouse gases on the European and regional (countrywise) scale. Two Lagrangian type of air pollution dispersion models have been adapted and were utilized in order to arrive at accurate description of process fluxes at the regional scale. The most important improvements in the ECN COMET model were: better and higher resolution 3D trajectory data and a better parametrisation of the mixing layer height. Furthermore measured vertical concentration gradients measured at the tall tower at Cabauw (NL) were used compared to the previous use of data for one measurement height only. The new version of COMET produces a much better comparison between concentration measurements and the model values, especially for the nighttime accumulation values. The linear regression of hourly data for the whole year of 2002 of methane between observations and the model delivers a correlation coefficient (r2) of 0.69. This means that the explained variation between measurements and model values has improved with almost a factor of two. From forward and inverse calculation we find that mixed layer averaged concentrations of methane at Cabauw are mainly determined by sources within The Netherlands. In the study period the contribution of The Netherlands to the methane in air at Cabauw in excess to the background is 65-80%. The total emission of the 6 areas of The Netherlands is estimated by the model inversions with a rather stable emission rate of 70-80% of the value of the estimate of the emission inventory METDAT, i.e. 1020 kton/yr for the base year 1995. The uncertainty analysis using a Monte Carlo error simulation technique indicates a rather high uncertainty of the individual emissions of the six separate emission areas in The Netherlands. The uncertainty is higher for the whole 2002 year period than for a shorter six week period in September/October 2002. After updating the emissions with the optimised values as calculated with the inverse COMET method, the correlation coefficient between observations and model values does not improve significantly. The Root Mean Square Error however is decreased substantially to about half the value obtained when using the uncorrected emissions. Using the optimized emission field also reduces the night-time overestimation of the modeled concentrations. The present results do not allow to discern the spatial distribution of the different source categories over The Netherlands. This can be attributed to several reasons like the error introduced by discretisation of the source areas, numerical uncertainties in the (big) solution space and the too high background concentrations used here (TM3, based on 1993 data).
Keywords Verification, Emissions, Greenhouse Gas, Methane, Nitrous Oxide, Transport Modelling, Lagrangian Dispersion Models, Regional Budget, Emission Inventory
2
ECN-C--04-027
INHOUD LIJST VAN TABELLEN
5
LIJST VAN FIGUREN
5
SYNOPSIS
6
SAMENVATTING
7
1.
INLEIDING
9
2.
DE TRANSPORTMODELLEN COMET EN FLEXPART 2.1.1 Voorwaartse berekening 2.1.2 Inverse berekening 2.1.3 Aanpassingen voorwaarts 2.1.4 Aanpassingen invers 2.2 FLEXPART 2.2.1 Modificaties
11 13 13 14 14 15 15
3.
RESULTATEN 3.1 Cabauw methaan concentratie metingen 3.2 COMET 3.2.1 Voorwaarts 3.2.2 Voorbeeld van de evaluatie van de regionale emissies voor Nederland met voorwaartse berekeningen 3.2.3 Invers 3.2.4 Resulterende emissieschattingen bij SVD inversie 3.3 Flexpart model resultaten
17 17 17 18
CONCLUSIES EN DISCUSSIE
30
4.
22 23 25 27
REFERENTIES
31
BIJLAGE A A.1
32 32 32 33 34 35 35 36 39 39 41 41 43
BIJDRAGE AAN JRC-ISPRA WORKSHOP OKTOBER 2003 The CHIOTTO project A.1.1 Scientific objectives and approach of CHIOTTO A.1.2 The CHIOTTO concentration monitoring program A.1.3 Rationale of the CHIOTTO measurement setup A.2 Cabauw long term tall tower observation of greenhouse gas vertical profiles A.2.1 Description and history of the ECN measurements at Cabauw A.2.2 Results A.3 Inverse modeling using tall tower observations A.3.1 Description of the COMET (inverse) transport model A.3.2 Description of the SVD method A.3.3 A demonstration of the method: COMET and CH4 References
BIJLAGE B
BIJDRAGE AAN JRC-ISPRA WORKSHOP OCT 2003
44
B.1
BOTTOM-UP INVENTORIES FOR USE IN INVERSE MODELS.
48
B.2
ATMOSPHERIC OBSERVATIONS
49
B.3
ATMOSPHERIC MODELS
51
B.4
LAGRANGIAN BACK TRAJECTORY MODELS
51
B.5
LAGRANGIAN PARTICLE DISPERSION MODELS
51
B.6
COUPLED MESOSCALE / GLOBAL MODELS
52
B.7
ZOOM MODELS
52
ECN-C--04-027
3
B.8
HIGH-RESOLUTION GLOBAL MODELS
53
B.9
INVERSION TECHNIQUES
53
4
ECN-C--04-027
LIJST VAN TABELLEN Tabel 3.1 Tabel 3.2 Tabel 3.3 Tabel 3.4 Tabel 3.5
Model berekeningen versus individuele meethoogten en gemeten menglaag concentratie ............................................................................................................ 20 Vergelijking van de fit tussen gemodelleerde en gemeten concentratie bij gebruik van verschillende emissie databases ......................................................... 21 Onderverdeling van Nederland in zes regio’s ........................................................ 23 Emissies in kton/jaar voor en na aanpassing met de TSVD inversie...................... 25 SVD berekende emissies van CH4 [kton/yr] voor enkele regios gebruikmakend van de Cabauw 2002 meetgegevens ....................................................................... 27
LIJST VAN FIGUREN Figuur 2.1 Schematische weergave van het Comet model ....................................................... 12 Figuur 3.1 De meetmast Cabauw ............................................................................................. 17 Figuur 3.2 Gemeten en gemodelleerde vertikale gradiënt van methaan (hoogtes 200, 120, 60, 20) in de mast van Cabauw. De zwarte lijn (3D_200) heeft de met COMET berekende concentratie weer, gebruikmakende van de 3D trajectorie informatie voor aankomsthoogte 200 meter. De gemeten concentraties worden weergeven met de andere lijnen (C_XXX, waarbij XXX staat voor de meethoogte in meters). ................................................................................................................... 18 Figuur 3.3 De modelberekening vergeleken met de gewogen concentratie op Cabauw.......... 19 Figuur 3.4 Correlatiecoefficient tussen berekende en gemodelleerde concentratie als functie van de uitsluiting van meetmomenten op grond van de standaarddeviatie in de concentraties van de onzekerheidstrajectorien [%] ...................................... 20 Figuur 3.5 Voorbeeld van de modelberekening waaraan een onzekerheid wordt toegekend aan de gemodelleerde concentraties op basis van de berekeningen op aankomstpunten rondom Cabauw. De grijze lijnen geven de modelmarges aan op basis van deze ‘onzekerheids’trajectoriën......................................................... 21 Figuur 3.6 Geoptimaliseerde model-concentratie.................................................................... 23 Figuur 3.7 Scatterplots van gemeten vs gemodelleerde concentraties in de periode sept-oct 2002. Links de directe voorwaarts berekende concentraties, rechts de door vermenigvuldiging van de SRM met de emissievector voorspelde concentratiereeks uitgezet tegen de gemeten waarden........................................... 24 Figuur 3.8 Concentraties in de menglaag, gemeten en gemodelleerd voor de periode 1 sep 2002 t/m 15 oktober 2002....................................................................................... 25 Figuur 3.9 Scatterplot van gemeten versus gemodelleerde concentraties met a posterior SVD emissies voor de zesweekse periode in 2002 .................................................. 26 Figuur 3.10 Scatterplot metingen vs COMET voorwaarts voor geheel 2002 ............................ 27 Figuur 3.11 Scatterplot metingen vs SRM invers voor SVD geheel 2002 .................................. 27 Figuur 3.12 Scatterplot metingen vs SRM invers na SVD voor geheel 2002 ............................. 27 Figuur 3.13 Met Flexpart gemodelleerde CH4 concentratie in Cabauw voor september 2002 vergeleken met de gemeten concentraties............................................................... 28 Figuur 3.14 Scatterplot van de gemeten methaan concentraties in Cabauw in september 2002 tegen de met Flextra gemodelleerde concentraties ....................................... 29
ECN-C--04-027
5
SYNOPSIS Cabauw CH4 CO2 COMET ECMWF EDGAR EMEP Flexpart Kyoto LPDM LOTOS N2O RMSE SF6 SRM SVD Trajectorie UN-FCCC
6
Meetlokatie meetmast KNMI, hoogte tot 216 m. boven NAP, nabij Schoonhoven Methaan Kooldioxide CO2-MEthane Transport model (ECN) European Center for Medium range Weather Forecasts te Reading, centraal Europees rekencentrum voor weersvoorspellingen RIVM emissie database voor broeikasgassen (Olivier & Berdowski, 2002) European Monitoring and Evaluation Program, vallend onder de Convention of Long-Range Transboundary Air Pollution (UN-ECE, 1979) Lagrangian Particle Dispersion Model van Universiteit München Plaats in Japan waar in 1997 klimaatverdrag werd gesloten Lagrangian Particle Dispersion Model Euleriaans verspreidingsmodel van TNO Lachgas Root Mean Square Error, maat voor de gemiddelde afwijking tussen schatting en werkelijkheid (hier tussen gemodelleerde en gemeten waarde) Zwavelhexafluoride Source Receptor Matrix Singular Value Decomposition, numeriekemethode voor matrixinversie Route die een luchtmassa aflegt United Nation Framework Convention on Climate Change
ECN-C--04-027
SAMENVATTING Doel van het project is de verbetering van de nauwkeurigheid van emissieverificatie op de lokale en Europees/regionale schaal. Twee Lagrangiaanse verspreidingsmodellen, COMET en FLEXPART, zijn aangepast en ingezet met als doel tot een nauwkeuriger procesbeschrijving te komen op de regionale schaal. De verbeteringen in het COMET model betreffen verbeterde trajectoriegegevens en een betere parametrisatie van de menglaaghoogte. Ook kan nu gebruik gemaakt worden van verticale concentratiegradiënten zoals gemeten in de Cabauw toren. Het COMET model presteert in de nieuwe versie veel beter in de vergelijking tussen voorwaarts berekende en gemeten concentraties. Vooral de accumulatie in de nachtelijke grenslaag wordt veel beter beschreven. Voor de uurlijkse gegevens over heel 2002 is de correlatiecoëfficiënt (r²) tussen meting en model voor methaan 0.69, dat wil zeggen dat zo’n 70% van de variatie in het signaal in de atmosfeer verklaard wordt door het model. Uit de voorwaartse en inversieberekeningen blijkt dat de voor de grenslaag gemiddelde concentratie van methaan in Cabauw voornamelijk bepaald wordt door bronnen in Nederland en wel voor 65-80%, afhankelijk van de periode. De totale emissie van de 6 onderscheiden brongebieden in Nederland is in alle inversieberekeningen redelijk stabiel en komt neer op een emissie van 70-80% (afhankelijk van de periode) van de METDAT emissiedatabase waarde van 1020 kton/jaar (basisjaar 1995). De Monte Carlo onzekerheidsanalyse van de SVD inversie geeft echter een vrij hoge onzekerheid aan van de berekende emissies van de 6 individuele Nederlands brongebieden. Deze onzekerheid is groter voor de periode van een vol jaar dan voor de periode van 6 weken in 2002. Na update van de emissies met door de TSVD berekende waarden neemt de correlatie tussen gemeten en gemodelleerde concentraties slechts weinig toe, de RMSE daarentegen wordt vrijwel gehalveerd, tot een waarde die ongeveer gelijk is aan de meetfout. De overschatting van de nachtelijke pieken wordt duidelijk teruggebracht na de update van de emissiewaarden. Het is moeilijk om met de huidige resultaten een beoordeling te geven van de ruimtelijke verdeling van de bronnen over Nederland. Redenen hiervoor zijn de (overigens voor de methode noodzakelijke) discretisering van de brongebieden en de onzekere numerieke oplossingsruimte door toepassing van de vijf broncategoriën en de enorme hoeveelheid concentratie metingen (10000 uurwaarnemingen in een jaar), die in combinatie met de grote spreiding in de orde van grootte van de contributiefactoren in de SRM leiden tot numerieke problemen in de matrixinversie. Mogelijk spelen ook de (te) hoge achtergrondsconcentraties (TM3, op basis van 1993emissies) in brongebieden een rol.
ECN-C--04-027
7
8
ECN-C--04-027
1.
INLEIDING
De emissies van de zes ‘Kyoto’ gassen (koolstofdioxide CO2, methaan CH4, lachgas N2O, fluorkoolwaterstoffen HFK’s, zwavelhexafluoride SF6) worden door de landen jaarlijks gerapporteerd volgens de door de UN-FCCC gegeven richtlijnen. In de schattingen van de emissies van de belangrijkste niet-CO2 broeikasgassen bestaan nog steeds grote onzekerheden. Voor methaan en lachgas wordt de onzekerheid in de totale emissie in Nederland geschat op 30-50% (Olivier et al., 2001). Wereldwijd zijn de niet-CO2 broeikasgassen verantwoordelijk voor 36% van het broeikaseffect; voor de Nederlandse emissies bedraagt dit (door de grote mate van gebruik van aardgas, een CO2-emissie arme brandstof) iets minder, namelijk 25%. Een aanzienlijk deel van de emissiereductiedoelstellingen op klimaatgebied is dus te behalen door het reduceren van de emissies van deze gassen. Vaak kunnen reducties van emissies van niet-CO2 broeikasgassen plaats vinden tegen lagere kosten per vermeden equivalentton CO2 dan bij het vermijden van CO2 emissies zelf. Bovendien kunnen de emissiereducties van de niet-CO2 gassen reeds op korte termijn effectief worden, hetgeen enige ruimte laat voor meer structurele en onvermijdelijke, maar meer ingrijpende verlagingen van de CO2 emissies op langere termijn. Om de onzekerheid in de emissieschattingen te verlagen en het monitoren van de emissieverlagingen van niet-CO2 broeikasgassen in de toekomst mogelijk te maken, is een onafhankelijke verificatie van die emissies noodzakelijk. Verificatie is een methode waarbij buitenluchtmetingen in combinatie met verspreidingsmodellen gebruikt worden om emissies te schatten. Dit is een door UN-FCCC erkende methode om, onafhankelijk van de traditionele methoden zoals gebruikt in de Emissieregistratie, te komen tot emissieschattingen. De nauwkeurigheid van verificatie wordt grotendeels bepaald door het vermogen van de verspreidingsmodellen om de metingen te simuleren (bijv: Builtjes et al., 2000). In een dichtbevolkt gebied als Nederland is het essentieel dat niet alleen de bijdragen van bronnen in Europa en de rest in de wereld goed beschreven worden, maar dat ook de bijdragen van lokale bronnen goed in kaart gebracht worden. In dit project worden twee lijnen van onderzoek gevolgd om de mogelijkheden tot verificatie van niet-CO2 broeikasgassen voor Nederland te verbeteren: De eerste lijn is het ontwikkelen van een variant van LOTOS voor de Nederlandse schaal, genaamd mini-LOTOS. De resultaten van dit onderzoek worden gepresenteerd in een parallel rapport, Roemer, 2004. Door de grotere resolutie in horizontale en verticale richting zal mini-LOTOS veel beter in staat zijn de waargenomen variaties in de gemeten concentraties van de niet-CO2 broeikasgassen in Nederland te simuleren. Dientengevolge zal verificatie van Nederlandse bronnen met een grotere betrouwbaarheid geschieden. De tweede lijn is een verbetering van de nauwkeurigheid op de Europees/regionale schaal. Twee Lagrangiaanse verspreidingsmodellen, COMET en FLEXPART, zullen worden aangepast en ingezet met als doel tot een nauwkeuriger procesbeschrijving te komen op deze schaal. De resultaten van dit onderzoek zijn gepresenteerd aan een groep van Nederlandse deskundigen op dit terrein. De belangrijkste bevindingen, aanbevelingen en opmerkingen van deze klankbordgroep zijn weergegeven in een derde rapport, Roemer en Vermeulen, 2004. Het belang van dit project is dat hieruit een oordeel ontstaat over de mogelijkheden van en de vereisten voor verificatie. Dit oordeel zal een basis vormen voor de inrichting van meetstations voor verificatie van broeikasgassen. Tevens kunnen de resultaten van dit project gebruikt wor-
ECN-C--04-027
9
den in de internationale discussie over verificatie. Deze discussie zal zeer gebaat zijn bij kwantificering van de mogelijkheden van emissieverificatie. De werkzaamheden van dit project komen neer op het verbeteren en toepassen van verspreidingsmodellen op de lokale schaal en regionale schaal om de invloed van Nederlandse bronnen voor de niet-CO2 broeikasgassen beter te simuleren, met als oogmerk daarmee (via dataassimilatie) de emissies van deze stoffen te verifiëren. Er wordt daarbij op twee schalen gewerkt: •
Lokale schaal: Het bestaande LOTOS model en de bijbehorende gegevensbestanden zullen worden aangepast naar domein en resolutie van het gebied; de beschrijving van de verticale uitwisseling zal worden herzien. Het aangepaste model (nu mini-LOTOS geheten) zal worden vergeleken met waarnemingen in Nederland (voor CH4). Voor Nederland zullen met behulp van dataassimilatie emissieschattingen worden gemaakt en uitgesplitst naar regio’s en brongroepen.
•
Regionale schaal: Het bestaande COMET model met bijbehorende gegevens zal worden verbeterd wat betreft meteorologische informatie (trajectoriën en menglaaghoogten) en de tijdsvariatie in de emissies. Voor de inversiemethode leidend tot emissieschattingen uit metingen, zal een geavanceerd Particle Dispersion Model (FLEXPART) geschikt gemaakt worden.
Op dit moment wordt op een aantal plaatsen in Nederland (Cabauw en Kollumerwaard) op continue basis methaan (en lachgas) gemeten. Om de regionale component in de modellen beter te kunnen toetsen zijn metingen op meer plaatsen nodig. De uitkomsten van dit project zullen van belang zijn bij een eventuele inrichting van een (klein) meetnet voor broeikasgassen. In dit rapport ligt de nadruk op de beschrijving van de werkzaamheden en resultaten van het werk aan het COMET en het FLEXPART model.
10
ECN-C--04-027
2.
DE TRANSPORTMODELLEN COMET EN FLEXPART
Het ECN COMET model (COMET: CO2 and Methaan Transport model) is een trajectoriemodel, dat gevoed wordt met meteorologische gegevens. Het model domein is Noord-West Europa. Tot nu toe zijn voor de modelberekeningen met COMET steeds zogenaamde terugwaartse 2D trajectorien gebruikt, berekend door andere instituten. Voor een hogere nauwkeurigheid zijn 3-dimensionale trajectorien echter veel beter. Met behulp van het programma FLEXTRA kunnen deze trajectorien nu door ECN zelf berekend worden voor ieder gewenst aankomst punt en lengte op basis van windvelden en andere meteoinformatie van het ECMWF te Reading. Door nu trajectorien te bepalen voor een aantal punten rond ieder aankomstpunt kunnen weerssituaties die leiden tot een grote onzekerheid in de trajectorieposities gedetecteerd worden en in de analyse uitgefilterd worden. De menglaaghoogte is een andere zeer belangrijke parameter in COMET. Bij analyses van de tot nu toe gebruikte menglaaghoogte gegevens bleek dat deze een grote overschatting maken. Door gebruik van meer geavanceerde methoden om de menglaaghoogten te berekenen uit meteorologische informatie (ECMWF), zoals de methode die gebruik maakt van een kritisch Richardson getal, is hier een grote verbetering mogelijk. Door verbetering van de beschrijving van de verandering van de emissies van methaan als functie van de tijd (seizoen, uur van de dag) in COMET, is een kleine verbetering van de voorspelde concentraties te verwachten. Het COMET model (CO2 MEthaan Transport) is ontwikkeld voor voorwaartse en inverse modelleringen. Het model gebruikt geen chemische module omdat dat niet nodig is voor methaan vanwege de lange verblijftijd in de atmosfeer en de tijdschaal van het transportmodel. COMET maakt gebruik van inverse trajectoriën. Eerdere berekeningen werden uitgevoerd door middel van trajectorie- en menglaag hoogtegegevens van EMEP MSC-W uit Oslo. De 96 uur inversie trajectoriën werden berekend door analyse van windvelden voor een hoogte van 975 mbar. Tijdens de ontwikkeling van het model bleek dat de accuraatheid van de voorspellingen uit het model voornamelijk bepaald werden door de emissies. Maar daarnaast is er een afhankelijkheid van: • de geografische resolutie van de methaanbronnen, • de veronderstelde achtergrondconcentraties, • de nauwkeurigheid van de meteorologische parameters, • de hoogte van de menglaag en het pad van de trajectorie. De maximale resolutie van het model is 1 boogminuut in breedtegraad en 1 boogminuut in lengtegraad. Dit correspondeert op de breedtegraad van Nederland met een celgrootte van het grid van 3 x 2 km. Dit maakt een zeer gedetailleerde beschrijving van de emissies mogelijk. Het COMET model onderscheidt twee verticale lagen, een menglaag en een reservoirlaag. Methaan emissies komen voor in de menglaag. De hoogte van de menglaag hangt af van de meteorologische condities en het seizoen. Deze kunnen worden ontleend aan directe metingen met radiosonde data of radar metingen of ze kunnen worden berekend met data uit een weermodel. In de berekeningen in dit rapport worden gegevens gebruikt van EMEP. Deze zijn ontleend aan een combinatie van weermodel gegevens en ballon observaties. De menglaag wisselt lucht uit met de reservoirlaag in een dagelijkse cyclus, onder invloed van het inkomende zonlicht. Als de hoogte van de menglaag verandert tijdens de dagelijkse cyclus, verdwijnt er lucht uit de menglaag naar of komt binnen uit de reservoirlaag. Hierdoor verdwijnt
ECN-C--04-027
11
een deel van het methaan uit de menglaag naar de reservoirlaag en vice versa. De berekeningen met het COMET model geven aan dat het methaan concentraties boven de achtergrondconcentratie, zoals die worden gemeten in Cabauw, vooral bepaald worden door emissies naar de menglaag tijdens de laatste 24-48 uur voor monstername. Dit wordt veroorzaakt door de dagelijkse verdunning van de menglaag met achtergrondlucht vanuit de reservoirlaag. Uitwisseling tussen de vrije troposfeer en de reservoirlaag is niet meegenomen, omdat de verwachting is dat dit een verwaarloosbaar effect heeft op de concentraties in de menglaag. De schatting van de methaanconcentratie bij de start van elke trajectorie is verkregen uit de wekelijkse gemiddelden van de berekende concentraties van het TM3 model van mei 1992– april 1993. Het TM3 model is een globaal drie dimensionaal Euleriaans model van de atmosfeer met een resolutie van 5.5*3 graden (Houweling et al., 1999). Langs de trajectorie worden alle emissies berekend die binnen de cirkel van een specifieke diameter vallen. De diameter van deze cirkel neemt lineair af van 150 km bij de start van het 96 uur trajectorie tot 20 km op de meetlocatie in Cabauw. Na berekening van de opnieuw berekende methaan concentraties wordt een nieuwe hoogte van de menglaag berekend. Het verlies of de winst van de opslaglaag is meegenomen en het model verplaatst het centrum van het ronde gebied met een nieuwe afgenomen diameter naar de volgende locatie op de trajectorie. Op de nieuwe positie wordt dezelfde berekening uitgevoerd door gebruik te maken van de daarvoor berekende concentratie. Deze procedure wordt herhaald tot het eindpunt wordt bereikt. Figuur 3.1 toont een schematisch overzicht van het model. Een meer gedetailleerde beschrijving van het COMET model is opgenomen in Vermeulen et al. (1997).
12:00
9:00
15:00 18:00
21:00
H
AH.98
Figuur 2.1 Schematische weergave van het Comet model In Figuur 3.1 is de route van een luchtpakket dat op 21:00 aankomt bij de meetmast in Cabauw weergegeven. Uitwisseling treedt op met een gebied in de cirkel onder de luchtmassa. De menglaaghoogte varieert met de tijd. De trajectoriegegevens kunnen op twee wijzen doorgerekend worden in het transportmodel COMET, in de zogenaamde voorwaartse en in inverse modus, beide rekenwijzen zullen in de nu volgende paragrafen toegelicht worden.
12
ECN-C--04-027
2.1.1 Voorwaartse berekening In de voorwaartse mode gebruikt het COMET model per stap in de tijd de emissies voor de grid cellen onder het huidige ronde brongebied. Het model berekent de verandering in concentratie en eventueel isotopen verhouding voor de gemodelleerde componenten in de luchtkolom. De hoogte van de kolom komt overeen met de menglaag. De emissie data komen voort uit de emissie inventarisatie per bron categorie op een regelmatig grid (Olivier et al. 1994; 1996).
2.1.2 Inverse berekening Inverse methoden worden voor verschillende componenten op mondiale schaal veel gebruikt om de bijdrage van verschillende bronsoorten of van verschillende regio’s te bepalen, gebruik makend van een wereldwijd netwerk van concentratiemetingen en een wereldwijd circulatiemodel (e.g. Enting et al., 1993). In onze studie gebruiken we een vergelijkbare methode op een kleinere ruimtelijke schaal waarin de emissiesterkte voor brongebieden rond Cabauw als parameters worden gebruikt. Metingen zijn beschreven als een functie van de parameters die moeten worden bepaald door het model. Normaal gesproken worden technieken toegepast zoals Kalman filteren (e.g. Zhang, 1997) om de voorspellingen van het model beter aan te laten sluiten op de observaties door middel van een a-priori set met gebiedsemissies. In onze berekeningen is geen a-priori informatie over de onbekende emissies gegeven en de set van vergelijkingen is direct opgelost. Het transportmodel is gebruikt om te bepalen wat de bijdrage van elk gebied is aan de gemeten concentratie. Met behulp van het model wordt voor ieder te onderscheiden brongebied en ieder tijdstip de bijdrage aan de concentraties in het aankomstpunt bepaald, dit levert een zogenaamde Source Receptor Matrix (SRM) op. Door deze SRM te vermenigvuldigen met de emisisie per brongebied (een vector) ontstaat een voorwaarste berekende voorspelde concentratie reeks in het aankomstpunt. Omgekeerd kan uitgaande van een gemeten concentratiereeks uit de SRM de beste schatting voor de afzonderlijke emissies bepaald worden. Dit probleem wordt hier opgelost door inversie van de SRM met behulp van een numerieke methode, de Singular Value Decomposition (SVD) (Press et al., 1992). Deze ontbinding van de SRM maakt het mogelijk om the rekenkundig slecht gedefinieerde componenten te identificeren en deze termen te verwijderen. De inversie routine is eerder gebruikt in het COMET model om 14C emissies van nucleaire reactoren van 14CH4 metingen in Cabauw af te leiden (Eisma et al., 1994, 1995). Nadere details van de methode worden toegelicht in Bijlage 1 van dit rapport. Om de fout in de berekeningen van het inversie model te schatten is een statistische Monte Carlo routine toegepast. In deze methode worden de elementen van de SRM verstoord met een willekeurige factor van 20-50%, in overeenstemming met een geschatte geaccumuleerde fout in de model SRM. Deze geaccumuleerde fout vertegenwoordigt de gecombineerde effecten van alle mogelijke foutenbronnen die eerder zijn genoemd. Op iedere verstoorde matrix wordt wederom de matrix inversie toegepast. Dit is verschillende keren herhaald met nieuw gegenereerde willekeurige fouten van de matrix elementen. De standaard deviatie in de berekende emissies per gebied is een maat voor de onzekerheid van het resultaat. Om tot betrouwbare emissieschattingen te komen moeten het aantal concentratiemetingen groter zijn dan het aantal gebieden waarvoor emissies worden berekend. In het inversiemodel zijn de emissies berekend voor gezamenlijke gebieden met het meer grove LOTOS ½°*1° grid (grofweg gelijk aan een 60 x 60 km2 grid). Nabij Cabauw wordt een fijne resolutie genomen te beginnen met het combineren van enkele cellen en verder bij Cabauw vandaan bevatten de combinaties meerdere cellen. De gecombineerde brongebieden lijken ook in vorm op de verwachte emissie of het landspatroon. De inversie methode staat toe dat de oplossing van de inversie wordt aangepast voor een gecombineerde bron of om verschillende gewichten te geven aan de belangrijkheid van de oplos-
ECN-C--04-027
13
sing voor verschillende brongebieden. Al dan niet gedeeltelijk fixeren van delen van de emissiekaart kan bijvoorbeeld worden gebruikt om verschillende gewichten te geven aan data, afhankelijk van de onzekerheden. In onze berekeningen werden de emissies van zee en oceaan gebieden en emissies van afgelegen gebieden vastgelegd op de huidige emissie inventarisatie gegevens. Dit verkleint het aantal gebieden waarvoor een oplossing moet worden gevonden en verkleint daarmee de onzekerheid in de daar bepaalde emissies.
2.1.3 Aanpassingen voorwaarts De meteorologische bestanden benodigd voor de berekening van de trajectoriën wordt verkregen door toepassing van een set scripts en Fortran programma’s die uitgevoerd worden op de computer van het ECMWF te Reading. Op het ECMWF worden alle modelgegevens bewaard in een reusachtig (magneetbanden)archief met een capaciteit van vele Terabytes, dit systeem heet MARS. Door middels van een door ECMWF beschikbaar gesteld gebruikersprogramma kunnen gegevens vanuit het archief opgehaald worden in de vorm van zogenaamde GRIB bestanden. De door ons aangepaste scripts en programmatuur zoals oorspronkelijk ontwikkeld door Wotowa (2004) genereren voor het door ons gekozen model-domein van NW-Europa drie-uurlijkse bestanden van ongeveer 2 Mbyte omvang met alle meteorologische informatie benodigd voor de modellen Flextra en Flexpart. Per jaar beslaat de omvang van de meteorologische gegevens ongeveer 6 Gbyte. Flextra wordt vervolgens gebruikt voor berekening van uurlijkse trajectoriën op basis van geanalyseerde 3D ECMWF meteorologie voor de gekozen aankomstpunten. Dezelfde meteorologische gegevens dienen als invoer van een eigen rekenprogramma dat voor het gehele domein de menglaaghoogtes berekend volgens een zogenaamd critical Richardson number benadering. Deze berekende menglaaghoogtes worden vervolgens afhaneklijk van trajectoriepositie en tijdstip toegekend aan de Flextra trajectoriegegevens. De emissies van methaan worden in de voorwaartse modus berekend met behulp van de METDAT emissiedatabase (Berdowski et al, 2000). De resolutie van de emissies in de METDAT database is 10x10 boogminuten. Er wordt onderscheid gemaakt naar 72 brontypes. Bij deze METDAT database is de variatie van de brontermen gegevens volgens tijdblokfuncties per categorie. De tijdblokfuncties zijn gegeven per maand van het jaar, dag van de week en uur van de dag.
2.1.4 Aanpassingen invers De resolutie van de SRM matrix was in de oude versie van COMET 1x1 boogminuut. Deze resolutie is te groot om voor een lange tijdreeks van een of meer jaren, uur bij uur voor het hele domein alle gegevens op te slaan. Daarom is de resolutie terguggebracht tot 10x10 boogminuten, een reductie dus met een factor 100, maar gelijk aan de resolutie van de METDAT emissiedatabase. Voor vergelijkbaarheid met de resultaten van de inversies zoals uitgevoerd in dit project door TNO met het LOTOS model is gekozen voor het doorrekenen van een SRM voor dezelfde indeling in landen en categoriën, dat wil zeggen 12 landen (waarbij Nederland onderverdeeld wordt in 6 subdelen [NW,N,NO,MW,M,MO,ZW,Z,ZO)) en vijf broncategoriën (energie, landbouw, natuur, landfills en industrieel). Overigens is een voorwaartse run vrijwel identiek aan een run voor inverse berekening. In het laatste geval wordt alleen extra per broncategorie en regio per trajectorie de bijdrage aan de concentraties in het aankomstpunt bijgehouden.
14
ECN-C--04-027
2.2
FLEXPART
FLEXPART is een zogenaamd Lagrangiaans Particle Dispersion Model (LPDM). Het maakt gebruik van dezelfde meteorologische informatie als het reeds genoemde FLEXTRA model, dat de trajectorieinformatie levert voor COMET. In FLEXPART wordt een serie ‘deeltjes’ aankomend in de meetpunten terugwaarts in de tijd gevolgd. Het gaat hierbij om een instelbaar aantal, meestal enkele duizenden, luchtdeeltjes. Het exacte pad van ieder deeltje hangt af van aselect bepaalde afwijkingen die samenhangen met de atmosferische verspreiding. Deze benadering is veel rekenintensiever dan de simpele benadering in COMET, waarbij langs 1 trajectorie steeds alle emissies worden genomen uit cirkelvormige gebieden rond de trajectorderposities. De diameters van deze gebieden zijn afhankelijk van de reistijd en hebben een van tevoren vastgestelde waarde. In FLEXPART hangt het feit of een bepaald gebied bijdraagt aan de concentratie in het aankomstpunt af van de daadwerkelijke actuele mate van atmosferische menging op ieder punt en tijdstip. De uitvoer van FLEXPART is direct toepasbaar voor inversiemethoden. Bij deze inversies worden uit de concentratiemetingen in de receptorpunten direct de verantwoordelijke emissiesterkten bepaald. Tevens leveren de inversies een (statistisch) getal voor de betrouwbaarheid van de bepaalde emissies op. Een belangrijk voordeel van een Lagrangiaans model is dat er geen kunstmatige numerieke diffusie optreed zoals dat in Euleriaanse modellen altijd het geval is. Dit is met name dicht bij een bron van belang of als de concentratie van een grondgebonden receptor gemodelleerd wordt. In een Euleriaans model wordt de emissie daar direct in de hele box van de actieve gridcel gemixt. De beschrijving van de turbulentie in FLEXPART lijkt op die gebruikt bij zgn Large Eddy Simulation en is meer accuraat dan in Euleriaanse modellen haalbaar is. Een deeltjes model heeft in principe een oneindig hoge resolutie. In de praktijk wordt deze echter beperkt door het aantal deeltjes dat losgelaten kan worden en door de resolutie van de meteorologische data. FLEXPART kan in voorwaartse mode worden gebruikt om de dispersie vanuit een bron te simuleren of in inverse mode om te berekenen welk gebied heeft bijgedragen aan een bepaalde gemeten concentratie. Voorwaartse en achterwaartse trajectorien kunnen worden berekend. In FLEXPART is een uitgebreid algoritme opgenomen waarmee het effect op de verspreiding door convectieve beweging in de atmosfeer kan worden meegenomen. Vanwege de lineaire opzet van het model is FLEXPART relatief gemakkelijk in te zetten voor inverse berekeningen. Ook kan voor elk deeltje in een luchtmassa worden gevolgd wat het brongebied is geweest. De toegevoegde statistiek over het aantal deeltjes in een uiteindelijk grid box maakt het evalueren van model fouten en het geven van een onzekerheid in de voorspelde concentratie niveaus mogelijk.
2.2.1 Modificaties Modelleren van de concentraties van broeikasgassen is niet standaard voorzien in Flexpart. De definitie van extra componenten is in principe echter vrij eenvoudig en ook voorzien in de modelopzet. Wel een probleem is dat Flexpart normaliter gebruikt wordt om voor slechts enkele bronnen de invloed op de concentraties te berekenen voor een korte tijdsperiode en met constante emissiesterktes. Bovendien geldt voor de meeste componenten dat hun verblijftijd in de atmosfeer slechts enkele dagen is, dit in tegenstelling tot de broeikasgassen, waar de verblijftijd in jaren tot honderden jaren geteld wordt. Voor het beschrijven van de verspreiding van emissies uit een domein dekkende emissie database is in Flexpart 5.0 voorzien door koppeling met de RIVM Edgar database (Olivier & Berdowski, 2001), dit is toegepast voor de berekening van NOx concentraties. Hierbij wordt de Edgar data ingelezen en per run omgezet in een RELEASES bestand dat per gridcel een bron definieert met een bepaalde emissiesterkte en –duur. Voor een lange run zou echter voor ieder uur een dergelijke set emissies per gridpunt gegenereerd moeten worden, hetgeen zou leiden tot
ECN-C--04-027
15
enorme RELEASES bestanden en een grote rekentijd en veel schijf- en leeswerk naar harddisk. Daarom is gekozen voor het genereren van deze emissiegegevens per tijdstap in het geheugen. Ingebouwd zijn de EDGAR V3.2 (Olivier & Berdowski, 2001) emissiedatabase voor CO2, CH4 en N2O. De resolutie van de EDGAR database is 1x1 graden, het aantal broncategoriën is 92 en het basisjaar is 1995. De uitvoer van Flexpart is verder uitgebreid met een binair GRADS type bestand dat de berekende concentraties over het modelgrid op 0.5x0.5 graad resolutie per drie uur op vijf modelhoogtes bevat, zodat eenvoudig met het GRADS pakket (http://grads.iges.org/grads/grads.html) alle resultaten uitgevoerd kunnen worden naar contourplots e.d.
16
ECN-C--04-027
3.
RESULTATEN
3.1
Cabauw methaan concentratie metingen
De metingen van CO2 en niet-CO2 broeikasgassen in de hoge mast van Cabauw (zie Figuur 3.1) bestrijken twee meetintervallen. In de periode van 1992-1997 is op de hoogste positie (200m) de concentratie van CO2, CH4 en N2O gemeten in half uurs intervallen. Nadat de Cabauw toren voor enkele jaren wegens groot onderhoud gesloten was, is een tweede meetreeks ingezet vanaf juni 2000 die liep tot mei 2004. In deze periode zijn halfuurlijkse vertikale gradiënten van CO2, CH4 en N2O gemeten op meethoogten 20, 60, 120 en 200 meter. De N2O metingen zijn halverwege 2002 wegens problemen met de precisie en uiteindelijk falen van de detector beëindigd. Volgens planning zal vanaf november 2004 een verbeterde meetreeks ingezet worden volgens de CHIOTTO (http://www.chiotto.org; zie Bijlage A) standaard met zeer hoge precisie.
EC Amsterdam Den
Utrecht
Cabauw Rotterdam
Figuur 3.1 De meetmast Cabauw
3.2
COMET
Door samenvoegen van de verschillende versies van COMET is het model systeem beter onderhoudbaar, omdat wijzigingen nu gelijktijdig plaatsvinden voor alle type berekeningen (voorwaarts, invers en trajectory coverage). Een bijkomend efficiëntie voordeel is dat ook voorwaartse en inverse berekeningen nu tegelijkertijd in dezelfde run kunnen plaatsvinden.
ECN-C--04-027
17
Voor het testen van de modelperformance is gezocht naar een periode van enkele weken die moest voldoen aan de volgende voorwaarden: • Grote dagelijkse variatie in de periode gekenmerkt door grote opbouw van concentraties in de nachtelijke stabiele grenslaag. • Grote synoptische variatie, gekenmerkt door opbouw van concentratie gedurende een periode van enkele dagen, meestal gevolgd door een snelle afbraak. • Hoge beschikbaarheid van metingen van de verticale concentratiegradiënt. In de periode 2000-2004 voldoet de tijdspanne van 1 september 2002 tot en met 15 oktober 2004 aan deze voorwaarden. Vanwege de lage precisie en kortere beschikbaarheid van N2O gradiënten is ervoor gekozen om alleen naar de in Cabauw gemeten verticale CH4 gradiënten te kijken. De gegevens van Kollumerwaard waren helaas nog niet beschikbaar voor de modeltest. Overigens is voor het testen van de modelperformance zeker interessant om te weten wat de mogelijkheden met slechts één meetpunt eigenlijk zijn.
3.2.1 Voorwaarts De gemeten gradiënt van methaan en de met COMET in voorwaartse modus voorspelde concentraties worden voor de genoemde periode getoond in Figuur 3.2. C_200 C_120
2.9
C_60
C_20 3D_200
2.7 2.5 2.3 2.1 1.9 1.7 21-9-02 0:00
26-9-02 0:00
1-10-02 0:00
6-10-02 0:00
11-10-02 0:00
Figuur 3.2 Gemeten en gemodelleerde vertikale gradiënt van methaan (hoogtes 200, 120, 60, 20) in de mast van Cabauw. De zwarte lijn (3D_200) heeft de met COMET berekende concentratie weer, gebruikmakende van de 3D trajectorie informatie voor aankomsthoogte 200 meter. De gemeten concentraties worden weergeven met de andere lijnen (C_XXX, waarbij XXX staat voor de meethoogte in meters). Duidelijk is dat de berekende concentratie grote gelijkenis vertoont met de gemeten concentraties, maar dat afhankelijk van de menglaaghoogte op ieder moment de gelijkenis met de concentratie op een of meer van de meethoogten varieert. In figuur 3.3 is daarom de gemiddelde concentratie van de atmosferische menglaag berekend door een gewogen gemiddelde te nemen van de gemeten concentraties tot de gemodelleerde menglaaghoogte in Cabauw. In Tabel 3.1 is het resultaat van lineaire regressie van voorspelde versus gemeten concentraties en de gewogen gemiddelde concentratie in de menglaag weergegeven. Duidelijk is dat gemiddeld genomen de 200 meter waarneming de laagste correlatie met de modelwaarden vertoont.
18
ECN-C--04-027
Van de waarnemingen geeft de concentratie op 20 meter hoogte de beste correlatie met de modelwaarden. Echter de beste correlatie wordt gevonden tussen de gewogen concentratie in de menglaag en de modelwaarden. Duidelijk uit deze cijfers is dat waarnemingen alleen op grondniveau of alleen op de top van een hoge mast voor dit modeltype en langlevende tracers als methaan niet voldoende zijn. Metingen van de verticale gradiënt op meerdere hoogten in de grenslaag blijken eigenlijk noodzakelijk.
Methaan concentratie in ppb
C MixL (metingen) 3.7
3D_200
3.2
2.7
2.2
1.7 21-9-02 0:00
26-9-02 0:00
1-10-02 0:00
6-10-02 0:00
11-10-02 0:00
Figuur 3.3 De modelberekening vergeleken met de gewogen concentratie op Cabauw Het COMET model veronderstelt dat de beweging van de hele menglaag en de reservoirlaag door één trajectorie beschreven kan worden. Deze versimpeling maakt het echter wel noodzakelijk dat de juiste keuzes gemaakt worden voor het type trajectorie dat gebruikt wordt om dit pad te beschrijven. De berekening van trajectoriën met Flextra biedt daartoe diverse opties. Een van de opties is het selecteren van een aankomsthoogte. Daartoe kan in het geval van Cabauw logischerwijze gekozen worden tussen 20, 60, 120 en 200 meter. In Tabel 3.1 (kolom 3D_200, 3D_120, 3D_60 en 3D_20) zijn voor de beschouwde testperiode de resultaten van COMET berekeningen voor 3D trajectoriën voor deze vier aankomsthoogten weergegeven. Het blijkt dat trajectoriën berekend voor een aankomsthoogte van 20 meter de beste resultaten leveren. Naast 3D trajectoriën, die op enig punt van de trajectoriën vanwege verticale bewegingen ook tot transport buiten de menglaag kunnen leiden, kan Flextra ook trajectoriën berekenen die geforceerd binnen de menglaag blijven, wederom te berekenen voor de vier aankomsthoogten. Ook de resultaten voor deze trajectoriën zijn tevens in Tabel 3.1 weergegeven. Ook voor menglaagtrajectoriën blijkt dat deze berekend voor een aankomsthoogte van 20 meter de beste resultaten leveren. De verschillen tussen de resultaten berekent met 3D en menglaag trajectoriën zijn slechts zeer gering.
ECN-C--04-027
19
Tabel 3.1 Model berekeningen versus individuele meethoogten en gemeten menglaag concentratie Model run Mix_200 Mix_120 Mix_60 Mix_20 Correlation with different heights r2 C200 0.04 0.03 0.04 0.05 r2 C120 0.13 0.12 0.14 0.18 r2 C60 0.42 0.44 0.50 0.51 r2 C20 0.46 0.49 0.54 0.57 Correlation with mixed layer concentration r2 0.57 0.60 0.66 0.69 Slope 0.45 0.43 0.46 0.53
3D_200 3D_120 3D_60
3D_20
0.04 0.13 0.37 0.39
0.03 0.12 0.38 0.42
0.03 0.12 0.42 0.47
0.05 0.15 0.52 0.56
0.49 0.47
0.52 0.43
0.57 0.42
0.68 0.48
MX_xxx = Menglaagtrajectorie met aankomsthoogte xxx meter gebruikt 3D_xxx = 3d trajectorie gebruikt met aankomst op Cabauw op xxx meter hoogte r2 = kwadraat van de correlatiecoefficient slope = helling van de fit tussen model en meting.
Om een indruk te krijgen van de onzekerheden in trajectorieposities door de onzekerheden in de gebruikte 3D windvelden kan Flextra ook zogenaamde onzekerheidstrajectorien bepalen. Hierbij worden op een op te geven aantal punten op een cirkel (diameter instelbaar) rond het aankomstpunt tevens trajectorien berekend. Op momenten dat de onzekerheden in de windvelden ertoe leiden dat de vijf trajectorien ieder een sterk verschillend afgelegd pad beschrijven in horizontale en of verticale richtingen, zal dit leiden tot meer of minder verschillende voorspelde concentraties. Als maat van de onzekerheid kan men de de standaardafwijking in de 5 voorspelde concentraties voor enig moment gebruiken. Vervolgens kan men dit als criterium gebruiken om de meting en modelvoorspelling voor dat moment wel of niet te gebruiken voor bijvoorbeeld een inverse berekening. 1 0.9
Correlatiecoefficient
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
10
20
30
40
50
60
StDev Conc [%]
Figuur 3.4 Correlatiecoefficient tussen berekende en gemodelleerde concentratie als functie van de uitsluiting van meetmomenten op grond van de standaarddeviatie in de concentraties van de onzekerheidstrajectorien [%]
20
ECN-C--04-027
70
In Figuur 3.4 is de relatie tussen de correlatiecoefficient tussen gemeten en berekende waarden weergegevens als functie van de onzekerheid in de concentratie van de onzekerheidstrajectorien. Uit de figuur blijkt dat de correlatie toeneemt bij toenemende onzekerheid. Een optimum wordt gevonden bij een relatieve standaardafwijking van 6-9%. Met dat criterium is de verklaarde variatie in het concentratiesignaal maximaal 92% (R2=0.92). Met dit criterium wordt wel 70% van de metingen verwijderd.
Methaan concentratie in ppb
Op basis van de verschillende modelberekeningen met de onzekerheidstrajectorien kan ook voor elke meting een onzekerheid in de berekende concentratie worden weergegeven. Een voorbeeld van de marge is te zien in Figuur 3.5.
MX_200
2.9
mod-std
C MixL (metingen)
2.7
mod+std 2.5 2.3 2.1 1.9 1.7 11-9-02 0:00
12-9-02 0:00
13-9-02 0:00
14-9-02 0:00
15-9-02 0:00
16-9-02 0:00
17-9-02 0:00
18-9-02 0:00
19-9-02 0:00
20-9-02 0:00
21-9-02 0:00
Figuur 3.5 Voorbeeld van de modelberekening waaraan een onzekerheid wordt toegekend aan de gemodelleerde concentraties op basis van de berekeningen op aankomstpunten rondom Cabauw. De grijze lijnen geven de modelmarges aan op basis van deze ‘onzekerheids’trajectoriën Voor de onzekerheid in de voorwaarts berekende trajectoriën is verder natuurlijk ook het gebruikte emissiebestand van belang. In de bovenstaande berekening is steeds de METDAT emissiedatabase (Berdowski, 2002) gebruikt, zonder variatie van de emissies in de tijd. In Tabel 3.2 worden de resultaten weergegeven als deze variatie wel doorgevoerd wordt. Het blijkt dat de correlatie tussen gemeten en gemodelleerde concentraties in dit geval zelfs enigszins afneemt, dit terwijl het verrekenen van de tijdsblokfuncties een omvangrijke overhead oplevert. Tevens is in Tabel 3.2 weergegeven wat de resultaten zijn als gebruik gemaakt wordt van een ietwat oudere emissiedatabase met veel lagere resolutie (1x0.5 graden), te weten de LOTOS database (Builtjes, 1990). De correlatie tussen gemeten en berekende concentraties is vergelijkbaar met die berekend bij gebruik van de METDAT database. Het lijkt erop dat bij deze combinatie van metingen en model de hogere resolutie van de emissiedatabase in het geval van METDAT geen voordeel oplevert. Tabel 3.2 Vergelijking van de fit tussen gemodelleerde en gemeten concentratie bij gebruik van verschillende emissie databases METDAT 0.68
METDAT Novar 0.69
LOTOS 0.68
METDAT excl. NL 0.29
Slope
0.47
0.44
0.52
0.61
RMSE (ppb)
6.1
6.9
5.5
6.2
r
2
ECN-C--04-027
21
3.2.2 Voorbeeld van de evaluatie van de regionale emissies voor Nederland met voorwaartse berekeningen In de voorwaartse berekeningen kunnen verschillende brongebieden aan en uit worden gezet. Door dit te doen kan de bijdrage aan het concentratiepatroon door het desbetreffende gebied eenvoudig uitgerekend worden. Door vervolgens de procedure uit te voeren waarbij een optimum wordt gezocht in de overeenkomst tussen de gemeten en berekende concentratie kan de emissie voor het gebied (hier Nederland) worden gevonden. Een voorbeeld van het resultaat van zo’n exercitie is weergegeven in Figuur 3.6. De blauwe lijn geeft gemeten concentratie op Cabauw (het gewogen gemiddelde van de metingen op de verschillende hoogten). De achtergrondconcentratie in de modelrun voor deze episode in 2002 lag rond de 1808 ppb. (Op basis van de TM3 model data voor 1993) De modelberekeningen laten zien dat 60 % van de toename van de CH4 concentratie ten opzichte van het achtergrondniveau afkomstig is van bronnen binnen Nederland. Dat is belangrijk omdat het aangeeft dat een belangrijk deel van het signaal dat we waarnemen inderdaad door Nederlandse emissies wordt bepaald. De overeenkomst tussen de gemodelleerde concentratie en de gemeten concentratie op Cabauw kan nu geoptimaliseerd worden door de emissies in Nederland met een factor te vermenigvuldigen. Bij die optimalisatie moeten we een keuze maken, namelijk bij welke menglaaghoogte we de modelberekeningen nog meenemen als we de Nederlandse emissie proberen te herleiden. De menglaaghoogte-routine geeft nu gegevens tot een minimale hoogte van 50 meter. Het optreden van een inversielaag aan de grond heeft een sterke invloed hebben op het concentratiepatroon in de onderste luchtlagen. Ondanks de aanzienlijke verbetering van de correlatie tussen de metingen en het model na implementatie van de nieuwe menglaaghoogte-module kan het model deze sterke inversielagen nog niet goed doorrekenen. Dit probleem kunnen we oplossen door de menglaaghoogte voor die gevallen uit de op Cabauw gemeten temperatuurmetingen te herleiden. Dit is nu nog niet geimplementeerd. De factor die we vinden hangt af van die keuze. Accepteren we alle menglaaghoogten boven de 100 m dan vinden we 0.65, boven 120m 0,7 en boven de 200m 0.83. Indien de in het model berekende menglaaghoogte onder de 100m komt dan kunnen we aan de temperatuurgegevens op Cabauw meestal duidelijk zien dat er een grondinversie aanwezig is. Deze is dan uit die metingen goed te bepalen. Dus dit lijkt voor de berekeningen nu een duidelijk minimum. Wordt de menglaaghoogte 200m of meer dan verandert de factor die we nodig hebben om de berekeningen te optimaliseren niet meer en blijft op 0.83 staan. De conclusie is dus dat we voor deze periode een emissie vinden equivalent met een jaarlijks niveau van 880 ± 100 miljoen kg CH4. Dat komt goed overeen met het verwachtte emissieniveau. Volgens de respectievelijke Milieubalansen bedroeg in 1995 de nationale emissie 1192 miljoen kg CH4, in 2002 899 miljoen kg, een reductie van 25 %. De factor die we hier vinden is gebaseerd op de hier gebruikte episode van 6 weken. Dezelfde soort analyse kan uitgevoerd worden voor een aantal andere episoden. De vraag is dan of we ook de jaar en seizoen variatie kunnen laten zien.
22
ECN-C--04-027
400
Methaan concentratie in ppb
4.5 4.0
C MixL (metingen)
Optimalisatie 70 %
MixH Accepted
Bijdrage NL %
350
50% lijn
3.5
300
3.0 250
2.5 200 2.0 150
100
1.0
50
0.5 0.0 31-8-02 0:00
NL contrinb %
Menglaag accept 1/0
1.5
0 5-9-02 0:00
10-9-02 0:00
15-9-02 0:00
20-9-02 0:00
25-9-02 0:00
30-9-02 0:00
5-10-02 0:00
10-10-02 0:00
15-10-02 0:00
Figuur 3.6 Geoptimaliseerde model-concentratie. De groene lijn in Figuur 3.6 geeft aan wanneer de menglaaghoogte bij aankomst voldoet aan het criterium (hier >120m). De oranje lijn geeft de bijdrage van bronnen binnen Nederland aan de toename in de CH4 concentratie weer. De rode en blauwe lijnen geven de gemeten concentratie en de berekende concentratie weer. Een optimale overeenkomst wordt verkregen bij een Nederlandse jaarlijks emissieniveau van 880 ± 100 miljoen kg CH4.
3.2.3 Invers Voor dezelfde periode van 1 september 2002 tot en met 15 oktober 2002 is tevens een inverse berekening uitgevoerd. Daartoe is gebruik gemaakt van de 3D trajectorie voor het aankomstpunt op 20 meter hoogte, aangezien deze de beste resultaten bij de voorwaartse simulatie liet zien. De emissies zijn ingedeeld in 25 brongebieden over Europa, waarbij de dichtstbijzijnde ‘buur’landen van Nederland (DK, UK, Ei, B, D, F, Lux, S, N) ieder een brongebied vormen. Nederland zelf is onderverdeeld in zes brongebieden, volgens de volgende indeling (zie Tabel 3.3). Tabel 3.3 Onderverdeling van Nederland in zes regio’s Naam Breedtegraad (NB) Lengtegraad (OL) NL-ZW 50,5 - 51,7 3,0 – 5,5 NL-ZO 50,5 - 51,7 5,5 – 7,5 NL-MW 51,7 – 52,5 3,0 – 5,5 NL-MO 51,7 – 52,5 5,5 – 7,5 NL-NW 52,5 – 53,5 3,0 – 5,5 NL-NO 52,5 – 53,5 5,5 – 7,5
Gebied Zeeland, Noord-Brabant west Noord-Brabant oost, Limburg Randstad Gelderland, Overijssel Noord-Holland Friesland, Drente, Groningen
Voor de emissies buiten deze gebieden levert de METDAT database de gegevens uit de EDGAR database. Deze emissies worden samengevoegd tot 9 gebieden (ZW, ZM, ZO, MW, MM, MO, NW, NM, NO) die samen het gehele domein van het COMET model beslaan. Voor ieder brongebied worden de emissies in vijf emissiecategorieën verdeeld: Industrie, Landbouw, Landfills, Natuur en Overig, zodat in totaal 125 brontypen ontstaan.
ECN-C--04-027
23
Het omzetten van de hoge resolutie METDAT database in een grovere regio emissie database leidt natuurlijk tot een iets andere voorspelde concentratiereeks. Door de met COMET berekende SRM te vermenigvuldigen met de gemiddelde emissies voor de 125 brontype en deze te vergelijken met de rechtstreeks berekende concentraties en de gemeten concentraties, kan beoordeeld worden hoe groot deze fout is. De correlatie tussen beide berekende concentratiereeksen is hoog (r2=90%), men verliest door deze aktie van het omzetten van de hoge resolutie METDAT database in 125 bronnen (25 gebieden en 5 emissiecategoriën) dus slechts 10% van de verklaarde variatie in het concentratiesignaal. De correlatie tussen de met de 125 brontypen berekende concentraties en de metingen is vrijwel gelijk aan die met de geheel METDAT gebruikende berekeningen (r2=0.63 resp. 0.602; N=1060). De RMSE tussen model en metingen neemt wel toe van 7.6 naar 9.1 ppb. De correlaties zoals hier gegeven zijn berekend door van alle concentraties de achtergrondwaarde af te trekken en bestaan dus uit de geschatte bijdrage van de bronnen exclusief de achtergrond. Deze correlatiecoëfficiënten vallen enigszins lager uit dan de eerder gepresenteerde getallen. In Figuur 3.7 worden de scatterplots van beide concentratie reeksen getoond uitgezet tegen de gemeten concentratiereeks. In Figuur 3.8 zijn de drie concentratiereeksen tezamen afgebeeld als functie van de tijd. Door de aggregatie naar grotere emissiegebieden neemt de breedte van de concentratiepieken toe. Door het uitsmeren van hoge emissies in de regio van de mast over een groter gebied neemt tevens de invloed van de lokale bronnen op het modelresultaat toe, waardoor de hoge pieken in stabiele situaties nog wat hoger berekend worden dan in de gewone voorwaartse berekening. . COMET forw ard performance R=0.794 R^2=0.630 Slope=1.513 Intercept=-0.04 RMSE=7.4 ppb
0
0
Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
1
1.1 1.2 1.3 1.4 1.5
1.1 1.2 1.3 1.4 1.5
SRM ante TSVD R=0.776 R^2=0.602 Slope=1.617 Intercept=-0.01 RMSE=9.1 ppb
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
1.1
1.2
1.3
1.4
1.5
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
1.1
1.2
1.3
1.4
Figuur 3.7 Scatterplots van gemeten vs gemodelleerde concentraties in de periode sept-okt 2002. Links de directe voorwaarts berekende concentraties, rechts de door vermenigvuldiging van de SRM met de emissievector voorspelde concentratiereeks uitgezet tegen de gemeten waarden.
24
ECN-C--04-027
1.5
Methaan 2
Methaan excess [ppm]
Measurement SRM conc Comet forw ard SRM post SVD
1
2-9-2002 3-9-2002 4-9-2002 5-9-2002 6-9-2002 7-9-2002 8-9-2002 9-9-2002 10-9-2002 11-9-2002 12-9-2002 13-9-2002 14-9-2002 15-9-2002 16-9-2002 17-9-2002 18-9-2002 19-9-2002 20-9-2002 21-9-2002 22-9-2002 23-9-2002 24-9-2002 25-9-2002 26-9-2002 27-9-2002 28-9-2002 29-9-2002 30-9-2002 1-10-2002 2-10-2002 3-10-2002 4-10-2002 5-10-2002 6-10-2002 7-10-2002 8-10-2002 9-10-2002 10-10-2002 11-10-2002 12-10-2002 13-10-2002 14-10-2002 15-10-2002
0
Figuur 3.8 Concentraties in de menglaag, gemeten en gemodelleerd voor de periode 1 sep 2002 t/m 15 oktober 2002
3.2.4 Resulterende emissieschattingen bij SVD inversie De afkapwaarde voor de eigenwaarde, een belangrijke factor bij de SVD matrixinversie techniek, is bij de SVD inversie vastgesteld op 10% van de maximum eigenwaarde. De emissies van de brongebieden Luxemburg, Frankrijk, Ierland en Zee zijn in deze versie gefixeerd op de METDAT waarde. In Tabel 3.4 zijn de ante en post SVD gemiddelde emissies van de 25 brongebieden weergegeven tezamen met de updatefactor. De emissie van Nederland en Duitsland wordt in de post SVD berekening met zo’n 20% verlaagd. Tabel 3.4 Emissies in kton/jaar voor en na aanpassing met de TSVD inversie Regio Emissie SVD Emissie Updatefactor % België 660 353 54 Duitsland 5369 4511 84 NL WM 187 112 60 NL ZW 161 42 26 NL ZO 174 172 99 NL MO 191 0 0 NL NW 112 33 30 NL NO 195 409 210 NL totaal 1020 769 74 De correlatie tussen berekende en gemeten concentraties is met de posterior emissie in redelijke mate verhoogd. De voorspelde concentratiereeks na toepassing van de emissies berekend met de SVD inversie zijn tevens weergegeven in Figuur 3.8. Duidelijk is te zien dat de match tussen metingen en berekeningen behoorlijk is toegenomen. De slope van de beste rechte is echter wel aanzienlijk van waarde veranderd, van 1.41 naar 0.82. De RMSE is gereduceerd naar 4.0 ppb. ECN-C--04-027
25
De resulterende vergelijking tussen gemeten en berekende concentraties is weergegeven in de scatterplot van Figuur 3.10.
0
Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1.1 1.2 1.3 1.4 1.5
SRMpost TSVD R=0.805 R^2=0.649 Slope=0.819 Intercept=-0.01 RMSE=4.0
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
1.1
1.2
1.3
1.4
1.5
Figuur 3.9 Scatterplot van gemeten versus gemodelleerde concentraties met a posterior SVD emissies voor de zesweekse periode in 2002 Na de succesvolle inversie van de zes weken in de testserie is tevens een inversie berekening uitgevoerd met COMET voor het gehele jaar 2002. In tabel 3.5 is het resultaat van deze berekening weergegeven. Ook in deze volle periode van een jaar wordt de CH4 emissie van Nederland door de SVD optimalisatie naar beneden bijgesteld, en wel met zo’n 10%. De correlatie tussen de voorwaartse berekening met COMET met de gemeten menglaagconcentratie is vergelijkbaar met die voor de zesweekse periode en dus ook vrij goed. De met de discrete brongebieden berekende concentraties geven in dit geval zelfs een hogere correlatie met de metingen. De optimalisatie met SVD levert voor het volle jaar een kleinere winst in de correlatiecoefficient op dan in het geval van de zesweekse periode (0.67->0.69). De winst in de RMSE in de inverse concentratiereeks is vrij groot (3.1->1.6). De scatterplots van de drie modelreeksen versus de meetresultaten zijn weergegeven in Figuur 3.10, 3.11 en 3.12.
26
ECN-C--04-027
SRM post TSVD R=0.833 R^2=0.693 Slope=0.82 Intercept=-0.005 RMSE=1.6 ppb
0
1 Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0
Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Modelled conc [ppm] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1.1 1.2 1.3 1.4 1.5
1.1 1.2 1.3 1.4 1.5
1.1 1.2 1.3 1.4 1.5
SRM ante TSVD R=0.821 R^2=0.675 Slope=1.50 Intercept=-0.018 RMSE=3.1 ppb
COMET forw ard performance R=0.816 R^2=0.666 Slope=1.18 Intercept=-0.001 RMSE=2.3 ppb
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
1.1
1.2
1.3
1.4
1.5
Figuur 3.10 Scatterplot metingen vs COMET voorwaarts voor geheel 2002
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
1.1
1.2
1.3
1.4
Figuur 3.11 Scatterplot metingen vs SRM invers voor SVD geheel 2002
1.5
0
0.1
0.2
0.3
0.4
0.5
0.6 0.7 0.8 0.9 1 Measured concentration [ppm]
Tabel 3.5 SVD berekende emissies van CH4 [kton/yr] voor enkele regios gebruikmakend van de Cabauw 2002 meetgegevens Regio Emissie SVD Emissie Updatefactor % België 660 1050 152 Duitsland 5369 4815 89 NL WM 187 133 71 NL ZW 161 128 79 NL ZO 174 159 91 NL MO 191 4 2 NL NW 112 184 164 NL NO 195 338 173 NL totaal 1020 948 92 Bij deze inversie is de afsnijwaarde van de eigenwaarden ingesteld op 1%. Dezelfde gebieden zijn vastgezet op de METDAT emissies als in de inversie van de 6-weekse reeks. De aangegeven trend in de updatefactoren is voor de meeste gebieden vergelijkbaar met die voor de 6weekse simulatie behalve voor België, waarvoor in de simulatie voor heel 2002 een hogere emissie aangegeven wordt. Opvallend is dat de emissie voor het brongebied NL Midden-Oost steeds naar een minimale waarde teruggebracht wordt. Gezien de verwachte variatie in emissies door het jaar heen zijn verschillen tussen de gemiddelde emissies over het hele jaar en die voor een 6-weeks periode te verwachten. Ook is de inversie voor de zes weekse periode duidelijk instabieler en is een hoge afkapwaarde noodzakelijk. De inversie voor een periode van een jaar is duidelijk stabieler dan inversies met het oude COMET model waar minder data (6-uurs ipv 1-uurs intervallen) en minder correlatie tussen metingen en model was (R2=0.4 i.p.v 0.7). Toevoeging van data van meer stations en een langere meetreeks zal ook hier het inversieresultaat nog stabieler maken.
Flexpart model resultaten
Gezien de veel grotere rekentijd is met Flexpart geen langere tijdreeks doorgerekend maar slechts een enkele maand, te weten september 2002. De rekentijd voor deze simulatie bedraagt enkele dagen op de beschikbare hardware (AMD Athlon 1 GHz). Het toepassen van de METDAT database in Flexpart heeft in dit project lange tijd tot grote praktische problemen geleid. Na succesvolle tests met de simulatie van 222Rn, CO2 en EDGAR CH4 bleken de simulaties ECN-C--04-027
1.2
1.3
1.4
Figuur 3.12 Scatterplot metingen vs SRM invers na SVD voor geheel 2002
In Tabel 3.5 zijn de met SVD berekende emissie voor enkele regios weergegeven ten opzichte van de uitgangsemissies.
3.3
1.1
27
1.5
met METDAT CH4 emissies steeds tot wonderlijke concentratiepatronen boven Europa te leiden. Uiteindelijk bleek dit terug te leiden op foutieve emissies voor bepaalde kolenmijnen die een factor 10000 groter zijn in de METDAT database dan andere bronnen. Na verwijdering van deze bronnen uit de database levert Flextra wel redelijke resultaten. Het probleem van Flextra met deze enorme bronnen komt voort uit het beperkte aantal deeltjes (2 miljoen) dat Flextra kan bewerken op onze hardware. Met een vaste massa per deeltje en enkele zeer grote bronnen leidt dat tot een voortdurende herverdeling van massa waarbij de invloed van de andere relatief kleine bronnen na verloop van tijd in een simulatie geheel verloren gaat en nog slechts de invloed van enkele grote bronnen gesimuleerd wordt, zodat meestentijds in Cabauw geen extra methaan berekend wordt en slechts af en toe enorme pluimen ‘overwaaien’. In Figuur 3.13 is afgebeeld de door Flextra voorspelde CH4 concentratie als functie van de tijd. Tevens is in deze figuur de gemiddelde menglaagconcentratie weergegeven. Het dynamische gedrag van de dagelijkse gang in de CH4 concentratie wordt door Flexpart redelijk beschreven in deze periode. Ook de synoptische variatie, de ‘onderkant’ van de dagelijkse pieken wordt redelijk goed gevolgd. De meeste nachtelijke pieken in concentratie worden door Flexpart onderschat. Dit geldt vooral in het midden van de beschouwde periode van 12 t/m 19 september. 2.7
2.6 Flextra CH4 Meas CH4 conc
2.5
2.4
CH4 [ppmv]
2.3
2.2
2.1
2
1.9
1.8
1.7 9-1-02 0:00
9-6-02 0:00
9-11-02 0:00
9-16-02 0:00
9-21-02 0:00
9-26-02 0:00
10-1-02 0:00
Date/time
Figuur 3.13
Met Flexpart gemodelleerde CH4 concentratie in Cabauw voor september 2002 vergeleken met de gemeten concentraties
De correlatieplot van de met Flexpart gemodelleerde concentraties versus de gemeten waarden is getoond in Figuur 3.14. De correlatiecoëfficiënt is aanzielijk lager dan bij de vergelijking van gemeten en de met COMET gemodelleerde waarden, nu is r2 namelijk 0.37 in plaats van ongeveer 0.70. Waarden als nu voor Flextra zijn overigens niet ongebruikelijk bij simulaties met modellen op regionale of continentale schaal, zoals LOTOS, TM5 e.d..
28
ECN-C--04-027
2.95
2.75
CH4 model [ppmv]
2.55
2.35 y = 0.4903x + 0.9891 2 R = 0.3683
2.15
1.95
1.75 1.75
1.95
2.15
2.35
2.55
2.75
2.95
CH4 measured [ppmv]
Figuur 3.14
Scatterplot van de gemeten methaan concentraties in Cabauw in september 2002 tegen de met Flextra gemodelleerde concentraties
ECN-C--04-027
29
4.
CONCLUSIES EN DISCUSSIE
De verbeteringen in het COMET model betreffen verbeterde trajectoriegegevens en een betere parametrisatie van de menglaaghoogte. Ook kan nu gebruik gemaakt worden van verticale concentratiegradiënten zoals gemeten in de Cabauw toren. Het COMET model presteert in de nieuwe versie veel beter in de vergelijking tussen voorwaarts berekende en gemeten concentraties. Vooral de accumulatie in de nachtelijke grenslaag wordt veel beter beschreven. Voor de uurlijkse gegevens over heel 2002 is de correlatiecoëfficiënt (r²) tussen meting en model voor methaan 0.69, dat wil zeggen dat zo’n 70% van de variatie in het signaal in de atmosfeer verklaard wordt door het model. Het blijkt hierbij onontbeerlijk om naast gegevens van de concentratie op grondniveau ook meetgegevens op grotere hoogte in de grenslaag te hebben. Uit de voorwaartse en inversieberekeningen blijkt dat de voor de grenslaag gemiddelde concentratie van methaan in Cabauw voornamelijk bepaald wordt door bronnen in Nederland en wel voor 65-80%, afhankelijk van de periode. De gekozen bronverdeling in 25 gebieden in Europa en vijf broncategoriën blijkt voor de inverse berekeningen geen gelukkige keuze. Van de 25 gebieden zijn er slechts enkele die een signifikante invloed op het signaal hebben en de onderverdeling in vijf categoriën die voor het belangrijkste brongebied, Nederland, niet in ruimte en tijd duidelijk gescheiden zijn, maken een stabiele inversie met de huidige gegevens moeizaam. De totale emissie van de 6 brongebieden in Nederland is in alle inversieberekeningen wel redelijk stabiel en komt neer op een emissie van 70-80% (afhankelijk van de periode) van de METDAT emissiedatabase waarde van 1020 kton/jaar (basisjaar 1995). De Monte Carlo onzekerheidsanalyse van de SVD inversie geeft echter een vrij hoge onzekerheid aan van de berekende emissies van de 6 individuele Nederlands brongebieden. Deze onzekerheid is groter voor de periode van een vol jaar dan voor de periode van 6 weken in 2002. Een synthese experiment met geconstrueerde concentratiegegevens laat zien dat de SVD methode correct werkt voor de set van gegevens, maar de Monte Carlo analyse laat zien dat de introductie van een kleine fout tot grote onzekerheden in de bepaalde emissies leidt bij de gekozen indeling van brongebieden en broncategoriën. Na update van de emissies met door de TSVD berekende waarden neemt de correlatie tussen gemeten en gemodelleerde concentraties slechts weinig toe, de RMSE daarentegen wordt vrijwel gehalveerd, tot een waarde die ongeveer gelijk is aan de meetfout. De overschatting van de nachtelijke pieken wordt duidelijk teruggebracht na de update van de emissiewaarden. Probleem bij de inversies is dat in enkele gevallen een niet realistische grote verschuiving van emissies van het ene naar het andere gebied binnen Nederland wel een verbetering van de match tussen model en meting teweegbrengt. Hoewel dit voor het totaal van de Nederlandse emissies niet van belang is maakt dit het moeilijk om met de huidige resultaten een beoordeling te geven van de ruimtelijke verdeling van de bronnen over Nederland. Redenen hiervoor zijn de (overigens voor de methode noodzakelijke) discretisering van de brongebieden en de onzekere numerieke oplossingsruimte door toepassing van de vijf broncategoriën en de enorme hoeveelheid concentratie metingen (10000 uurwaarnemingen in een jaar), die in combinatie met de grote spreiding in de orde van grootte van de contributiefactoren in de SRM leiden tot numerieke problemen in de matrixinversie. Mogelijk spelen ook de (te) hoge achtergrondsconcentraties (TM3, op basis van 1993-emissies) in brongebieden een rol.
30
ECN-C--04-027
REFERENTIES Berdowski, J.J.M., G.P.J. Draaijers, L.H.J.M. Janssen, J.C.Th Hollander, M. van Loon, M.G.M. Roemer, A.T. Vermeulen, M. Vosbeek, H. Visser (2001): Sources, regional scaling and validation of methane emissions from the Netherlands and NorthWest Europe. Bilthoven, RIVM, report 410200084. Eisma, R., A.T. Vermeulen, K. van den Borg (1995): Radiocarbon 37, 2, 475-483. Olivier, J.G.J. and J.J.M. Berdowski (2001): Global emissions sources and sinks. In: Berdowski, J., Guicherit, R. and B.J. Heij (eds.) "The Climate System", pp. 33-78. A.A. Balkema Publishers/Swets & Zeitlinger Publishers, Lisse, The Netherlands. ISBN 90 5809 255 0. Press, W.H., B.P. Flannery, S.A. Teukolsky en W.T. Wetterling (1989): Numerical Recipes in Pascal. Cambridge University Press. ISBN 0521387671. Roemer, M. (2004): Methode ontwikkeling voor verbetering van de verificatie van niet-CO2 broeikasgassen – lokale schaal -. Apeldoorn. TNO, rapport R2004/501. Roemer, M. en A.T. Vermeulen (2004): Modelverbetering ten behoeve van verificatie van nietCO2 broeikasgasemissies – Samenvatting en discussie – Apeldoorn. TNO-MEP, rapport R2004/511. Vermeulen, A.T., R. Eisma, A. Hensen, J. Slanina (1999): Transport model calculations of NWEuropean methane emissions. Env Sci & Policy 2, 315-324. Wotowa, G., H. Kröger (1999): Testing the ability of trajectory statistics to reproduce emission inventories of air pollutants in cases of negligible measurement and transport errors. Atm Env 33, 3037-3043.
ECN-C--04-027
31
BIJLAGE A BIJDRAGE AAN JRC-ISPRA WORKSHOP OKTOBER 2003 CHIOTTO - CONTINUOUS HIGH-PRECISION TALL TOWER OBSERVATIONS OF GREENHOUSE GASES
A description of the project, a preview of the data it will provide and the use of its data in inverse modelling A.T. Vermeulen*, A. Hensen, M. Gloor2, A. Manning2, P. Ciais3, R. Eisma, W.C.M. van den Bulk, J.J. Mols, J.W. Erisman ECN – Air Pollution & Climate Change group, P.O. Box 1, 1755 ZG Petten, NL 2 Max-Planck-Institut für Biogeochemie, Postfach 10 01 64, 07701 Jena, D 3 CNRS-LSCE, Centre d’études l’Orme des Merisiers, Bât. 709, 91191 Gif-sur-Yvette Cedex, F * e-mail:
[email protected] A.1
The CHIOTTO project
A.1.1 Scientific objectives and approach of CHIOTTO The CHIOTTO project (http://www.chiotto.org) objective is to build an improved infrastructure for the continuous monitoring of the concentrations of greenhouse gases on the European continent above the surface layer using tall towers. This will be an important step towards a fully operational continuous observing system in the framework of the Kyoto Protocol for the sources and sinks of the most important greenhouse gases (CO2, CH4, N2O, CO, SF6) over Europe. An important aspect of the objective is the establishment of high quality calibrations for the existing and new atmospheric measurement stations, and the implementation of a near-online datatransmission system for tall tower measurements.
Figure A4-1 The influence function for the year 2002 of the 8 CHIOTTO tall towers derived by the COMET trajectory model; The white stars indicate the tall tower positions. In CHIOTTO the intercomparability of the concentration measurements between the institutes operating the air sampling networks will be monitored. Quality controlled atmospheric concen-
32
ECN-C--04-027
tration, CO2 flux and additional meteorological data will be archived in a data center accessible to the scientific community through the World Wide Web. We will integrate existing flux towers in the vicinity of the tall towers with the atmospheric stations networks in a synergetic approach enabling the tall towers to become atmospheric monitoring sites for use in transport models.
A.1.2 The CHIOTTO concentration monitoring program Table A1 Precision and accuracy goals for CHIOTTO Gas species CO2 CH4 CO N2O SF6 O2/N2 Rn δ13C-CO2 δ18O-CO2
Intra-laboratory instrumental precision 0.05 ppm 2.0 ppb 1.0 ppb 0.1 ppb 0.1 ppt 5 per meg 0.2 Bq m-3 or 10% ** 0.007 ‰ 0.03 ‰
Inter-laboratory calibration scale accuracy 0.10 ppm * 3.0 ppb 3.0 ppb 0.2 ppb 0.2 ppt 10 per meg 10% ** 0.01 ‰ * 0.05 ‰ *
* These values are the WMO/GAW goals for global network accuracy among different laboratories in the northern hemisphere. (In the case of the other species, no official WMO/GAW goals exist at this time). ** For Radon, these values are the detection limit goal.
In Europe tall towers in the Netherlands, Hungary, Germany and Sweden already are equipped with devices that permit online CO2 concentrations and in some cases related tracers. In the CHIOTTO-project these towers are extended with concentration measurements of more greenhouse gases and with new towers in Great-Britain, France, Poland and Italy to increase the spatial coverage of the studied area. The influence regions of these towers for the year 2002 together cover Europe very well (see Figure A4-1). This influence region analysis was performed using the COMET v2 Lagrangian transport model as described in the second part of this chapter. The area in Figure A4-1 with an influence larger than 0.2 can be considered as the area from which the sources and sinks have a significant influence on the measured concentrations at the towers. This analysis shows that this area covers most of Western Europe and extends up to about 500-1000 km from each individual tower. Table A2 Overview of the main parameters of the CHIOTTO tall towers Name Cabauw Griffin Hegyhatsal Orleans Norunda Florence Ochsenkopf Bialystok
NL UK H F S I D PL
Height (m) 200 232 117 203 102 245 177 330
Position Lon Lat 04°56’
51°58’
-2°59’ 16o39’ 2°07’ 17°28’ 11°16’ 11°49’ 22°45’
56°33’ 46o57’ 46°58’ 60°05’ 43°49’ 50°03’ 52°15’
CO2 4 1 4 4 4 1 3 5
Concentration measurement (levels) CH4 N2O SF6 CO 222Rn Flasks 3 4 4 4 4 1 1 1 1 1 1 3 1 1 1 1 3 4 4 4 4 1 1 1 1 1 1 3 3 3 3 1 3 5 5 5 5
Flux meas CO2 CH4 2 2 2
1
Operator ECN UEDIN ELTE LSCE LUPG UNITUS MPIBGC MPIBGC
The concentration measurements at the tall towers will have to be very precise and accurate, in order for them to be used successfully in transport model inversion studies to derive emissions. The high demands imposed on the measurements are summarized in Table A1. The same precision and accuracy targets have been accepted for application in the whole CarboEurope Integrated Project (IP), that has recently started in the framework of the EU FP6 research program, and of which the CHIOTTO project will become an integral part after the finalization of the FP5 period. A quick overview of the main characteristics of the tall tower stations of the CHIOTTO project
ECN-C--04-027
33
is given in Table A2.
A.1.3 Rationale of the CHIOTTO measurement setup Without knowledge of today’s sources (and sinks) of the most important greenhouse gases, theoretical predictions of future levels cannot be verified. Furthermore, without a sound and credible scientific understanding, it is very difficult to gain consensus to move forward on policy. Continued uncertainty about for example the importance of carbon sinks or the exact national greenhouse gas emissions will further limit the implementation and use of instruments aimed at reducing these emissions. In recent years it has become clear that the traditional atmospheric greenhouse gas observing system has severe gaps. Knowledge of for example the global carbon derived from atmospheric data was based on sparse sampling at ‘remote’ stations located far from the continents. This marine network is good enough to infer CO2 fluxes at the hemispheric scale, but cannot resolve these fluxes on the longitudinal scale. An example of the controversy to which this can lead is the current debate on the apportionment of the Northern Hemispheric CO2 sink between Northern America and Eurasia. The (lower) troposphere acts as an integrator (mixer) of spatially and temporally varying surface fluxes. The distribution of CO2 in the atmosphere and its time evolution can thus be used to quantify surface fluxes. However, the atmospheric approach to derive fluxes from observations of concentration on a regional scale needs sites better adapted to the specific circumstances at the regional level, where the precision and resolution in time and space of the measurements should match that of the major processes involved. Regional carbon budgets are currently constructed from approximately 100 points worldwide. This renders retrievals of the spatio-temporal distribution of surface fluxes highly uncertain (Rayner et al., 1999; Fan et al, 1998; Gloor et al, 1999). In addition, the understanding of processes from local studies [e.g. Feely et al, 1999; Wofsy et al, 1993) is difficult to ‘connect’ to global CO2 patterns because the undersampling of the atmosphere (Braswell et al, 1997). Without comprehensive spatial coverage, uncertainty cannot even be localized to transport model error, data error, or inversion procedures unequivocally. Until recently, continental monitoring stations have not been considered relevant because of the large variability of the signals, caused by the proximity to the land-biosphere atmosphere exchange fluxes and the large, spatially concentrated sources. Measurements in the free troposphere are to characterise the fluxes at the continental scale (e.g. Eurasia vs. N. America) because they provide a constraint of the annual mean and seasonal distribution of atmospheric concentrations over the interior of the continents. The horizontal gradients of greenhouse gas concentrations, which carry the information on the magnitude and spatial distribution of sources and sinks, are quite small. If we are to infer fluxes at the regional level, it is thus necessary to sample close to the earth surface and on a continuous basis to capture the signal of greenhouse gas exchange fluxes. This calls for measurements in the boundary layer. Here the variability in concentrations (diurnal cycles) is huge, because the air is to a large extent influenced by local sources. In order to separate the effect of local (few tens of km) variability from the regional signal, one needs to continuously monitor concentrations above the surface layer (100 meters), complemented for CO2 by eddy flux towers to characterize the contribution of ‘local’ biospheric exchange. We can take advantage of the strong mixing in the lowermost 1 to 2 km of the atmosphere during daytime caused by cumulus convection. If the gases are measured at sufficient height above ground (ideally a few hundred meters), then a fairly homogeneous signal that integrates fluxes over a footprint on the order of a circle of 500 to 1000 km is obtained (Gloor et al, 2001). Furthermore continuous measurements permit to optimise the signal to ‘noise’ ratio of the measurements. As the variance of the measured concentrations is the most relevant quantity for the precision achievable with inverse calculations, the number measurements should be as high as possible. Because of the natural variability of the signals in the atmosphere, the time series needed for successful inverse calculation of fluxes by atmospheric transport models (ATM’s) should extend over a long time frame of at least several years. Time series of 5-10 years or longer will be 34
ECN-C--04-027
required. In order to be able to distinguish trends in time, because of a response to e.g. climate change on ecosystems CO2 fluxes or emission reduction measures on CH4 emissions or land use changes on N2O emissions. At this time no infrastructure yet exists that will provide us with that data. A relatively new approach is to sample continuously CO2 and other greenhouse gases (some also related to the carbon cycle) like CH4, CO, N2O and tracers that help validate the realism of transport simulation with atmospheric models like SF6, C2Cl4 and 222Rn on tall towers. An investigation of the signal to ‘noise’ in highly resolved models of atmospheric transport indicates that tall tower stations spaced at approximately 500 to 1000 km from one another are needed to obtain sufficient statistics for robust conclusions on carbon fluxes (see Figure 1)
Figure A4-2 Model simulated standard deviation in the daytime CO2 concentration, representing the size of daily CO2 signal due to biogenic activity, as a function of height and longitude (from: Gloor M, et al., 2001). In the future remote sensing techniques will enable concentration measurements for CO2, CO, CH4 and maybe also other greenhouse gases using satellite observations. Main advantage of these techniques will be that they will provide excellent spatial coverage. Main problems of the remote sensing data for use in (inverse) modelling of earth-atmosphere fluxes are the inherently limited accuracy of the measurements, the limited (e.g. cloud-free) conditions for which observations will work and the limited vertical resolution (especially in the troposphere) of the data, in many cases only column totals can be derived. Development of these techniques will also still take quite some years from now before becoming operational. When these techniques become available they will certainly complement the ground-based measurements. On the other hand the ground-based observations will then also be useful as ground-truthing data for calibration and verification of the remote-sensing data.
A.2
Cabauw long term tall tower observation of greenhouse gas vertical profiles
A.2.1 Description and history of the ECN measurements at Cabauw The ECN measurements take place on a 213-meter high meteorological tower near the village of Cabauw, located in the center of The Netherlands (51°58’16’N, 04°55’36’E). Figure A4-3 shows the location of the site. The tower is owned by the Royal Netherlands Meteorological Institute Meteorological (KNMI). The site is located in an agricultural region with mainly grasslands. At about 40 km to the west a variety of sources (traffic, waste water treatment plants, ECN-C--04-027
35
etc.) are located in a city belt that runs from the Rotterdam area to the north in Amsterdam. The analytical equipment for the greenhouse gas measurements is placed in the basement of the tower. Ambient air is sampled continuously from 200 m, 120 m, 60 m and 20 m height to the basement, by a membrane pump at approximately 6 l.min-1 through an inlet filter, Permapure dryers, PTFE tubing and a high-efficiency particulate air filter, respectively. The GC system consists of a Carlo Erba 8000 GC equipped with FID and ECD detectors for detection of CH4 and N2O respectively. The chromatograms are analysed using Shimadzu ClassVP version 4.3 software. Every five minutes a new sample is injected according to the following scheme: high standard, low standard, 200m, 120m, 60m, 20m. This scheme is repeated every half hour. The working standards are calibrated against NOAA certified standards regularly. CO2 is measured using a Siemens Ultramat 5 NDIR. The resolution of this monitor is 0.1 ppm. Every 24 hours the monitor is autocalibrated using zero air and a working standard. This working standard is calibrated against a NOAA certified standard regularly. Meteorological data are provided by the KNMI. With steps of 20 m wind speed, wind direction, pressure, temperature, radiation and relative humidity are available as 30 min averages. All measurements are transmitted to ECN on a half hourly basis. Until 1997 CO2, N2O and CH4 were sampled only from the 200m level (Vermeulen et al, 1997). After refurbishment of the tower, measurements restarted in the year 2000 by sampling at the four levels as described before. L o p ik P o ld er are a
G ra ssla n d
1 k ilo m e te r
O rch ard s a n d G ra sslan d
2 0 0 m to w e r KNMI R iv e r L e k
ECN A m ste rd a m
T he hague
U tre ch t
C abauw R o tterd a m
Figure A4-3: Cabauw tower location and photograph. The photo shows only 3 of the 20m spaced platforms of the tower. The total tower height is 213 meter.
A.2.2 Results Figure A4-4 shows the full dataset available for CO2 at the 200 m level from October 1992 until December 2002. This dataset clearly shows the continuous increase in the CO2 concentrations over the Netherlands.
36
ECN-C--04-027
CO2 92-93 CO2 95 CO2 97 CO2 00 CMDL Ocean station M
490
CO2 concentration in ppm
475 460
CO2 94 CO2 96 CO2 01 CO2 02 Izana I.N.M.Spain
445 430 415 400 385 370 355
01-04
01-03
01-02
01-01
01-00
01-99
01-98
01-97
01-96
01-95
01-94
01-93
01-92
340
Figure A4-4: The CO2 concentration at 200m height at Cabauw in the period 1992 to 2002, compared with the background stations Izana (Ternerife) and Ocean station M (Neas Shetlands) In Table A3 the monthly means for the measurements in 2001 for all three gases and measurement height are shown. The seasonal variation (Figure A4-5) of the CO2 concentration pattern shows minimum concentrations in summer when uptake by vegetation occurs and the mixing layer height is generally high (in the order of 1000 to 2000 m). In wintertime respiration from vegetation, and emissions from anthropogenic sources increase the concentration and this increase is further enhanced because the emissions take place in a smaller air volume, with mixing layer heights that are generally in the order of 200 to 1000 m. In addition to the Cabauw dataset shown in the graph with an 3 hour average resolution, the data obtained by CMDL at Ocean Station M is shown. The trend and the seasonal variation agree on average (source: GAW database, T.Conway, CMDL and the Norwegian Meteorological Institute). The Cabauw station shows enhanced concentrations compared to those observed at ocean station M in winter, being within the source area. In spring and autumn periods occur where the Cabauw concentration levels decrease below the ocean station M levels. The uptake of CO2 by biota and the increased mixing layer as the air passes over the land cause this. Table A3 Monthly averaged concentrations at the Cabauw tower at different heights in 2001 All concentrations in ppm Month CO2 CO2 20 m 60 m Jan 409.5 406.7 Feb 395.6 393.5 Mar 397.4 396.1 Apr 387.8 386.3 May 385.2 382.9 Jun 381.0 377.1 Jul 378.3 372.6 Aug 380.4 373.2 Sep 374.7 371.5 Oct 389.6 386.4 Nov 395.0 388.3 Dec 399.4 392.8 Nr of meas 15341 15341 % coverage 88% 88%
ECN-C--04-027
CO2 120 m 403.8 393.3 395.6 386.0 382.6 376.4 370.0 372.0 370.1 385.5 386.5 390.5 15342 88%
CO2 200 m 395.9 390.0 392.5 383.9 379.2 373.8 366.8 367.3 366.9 379.4 381.3 386.9 15342 88%
CH4 20 m 2.205 2.098 2.071 1.946 1.928 1.976 2.026 2.108 1.998 1.972 2.182 2.207 11473 65%
CH4 60 m 2.163 2.041 2.008 1.908 1.897 1.910 1.954 2.010 1.942 1.939 2.030 2.051 11468 65%
CH4 120 m 2.125 2.034 1.999 1.893 1.888 1.894 1.917 1.968 1.905 1.909 1.962 1.988 11473 65%
CH4 200 m 1.993 1.961 1.952 1.875 1.871 1.878 1.882 1.916 1.877 1.857 1.908 1.936 11665 67%
N2O 200 m 0.313 0.311 0.309 0.309 0.312 0.312 0.312 0.312 0.312 0.310 0.308 0.310 13345 76%
37
42 0 20 m 60 m 120 m 200 m
CO 2 Concentration (ppm)
41 0 40 0 39 0 38 0 37 0
31-Dec-02
1-Jul-02
31-Dec-01
1-Jul-01
31-Dec-00
1-Jul-00
1-Jan-00
36 0
Figure A4-5: Seasonal averaged concentrations of the CO2 vertical profiles at Cabauw tower.
Figure A4-6 Time series of the CH4 concentration at 200m height at Cabauw in the period 1992-2002. The background monitoring stations Ocean station M and Izaña clearly show the north-south gradient in the methane concentration. For methane the data set for 200m is shown in Figure A4-6. In this graph also the data obtained from the Global Atmospheric Watch database for Ocean station M and for Izaña area plotted. In contrast to CO2, we do not expect significant periods where the CH4 concentration is below the background levels. From the graph it is clear that sometimes the Cabauw data shows concentration levels below the values observed at Ocean Station M, but almost never concentrations that are below the Izana data. Further evaluation of this comparison requires the information of the origin of the airmass that is actually sampled at a given moment in time. A first rough hypothesis based on the graph would be that the 2000-2001 data seem to approach the low background concentrations more frequently than in the period 1992-1994. This might indicate that CH4 emission reduction measures, like those already implemented in landfill management, actually are detectable already at this moment. Another possibility is that the station standards or the GC performance have changed slightly. This is subject of ongoing evaluation of the data by recalibration of the currently used NOAA calibration cylinders.
Well mixed conditions 38
ECN-C--04-027
Both the gradient data obtained in the period 1992-1997 and the extended gradients in the period 2000-2003 show that for CO2 in general during daytime the concentrations at the low levels and the 200m levels are close together with differences in the order of a few ppm or below 1 ppm. This is a well mixed situation this data will be the main data to use for large scale modeling and inverse techniques. The data obtained at the 200 m level is representative for a large area and the effect of ground level sources will be distributed though the whole mixing layer. As an example, in Figure A4-7 most daytime situations in August 2000 show a well-mixed situation for the CO2 data.
Strong gradient conditions During most nights and during some days, mainly in winter, there is a strong gradient and concentrations at all heights are elevated. Detailed information about the mixing layer height is needed in order to evaluate these data. Since in general these situations occur with relatively low wind speeds, the effect nearby sources will have a significant contribution. Until now no proper modeling tool exists to assess these situations but since they occur relatively frequent the development of these tools will be part of new and running projects. In Figure A4-7 for example the data for August 19 shows a strong gradient that clearly extends to an altitude above 200m. C O 2 con ce n t r a t ion a t C a ba u w Au gu s t 2000 490.00 20 m et er 60 m et er
470.00
120 m et er
Concentration [ppm]
450.00
200 m et er
430.00 410.00 390.00 370.00 350.00 1 7 /0 8 /0 0
1 8 /0 8/0 0
1 9 /0 8 /0 0
2 0 /0 8 /0 0
2 1/0 8 /0 0
2 2 /0 8 /0 0
23 /0 8 /0 0
2 4 /0 8 /0 0
2 5 /0 8/0 0
2 6 /0 8 /0 0
2 7 /0 8 /0 0
Da te
Figure A4-7: Half hourly values of the vertical concentration gradient of CO2 at Cabauw tower for a period of 10 days in August 2000. The measured concentration gradient of CH4 shows nearly the same pattern
Low inversion events Furthermore there is the situation where the mixing layer has a temperature inversion below the 200m level. In that case a strong gradient occurs at the low levels but the 200m concentration remains relatively low. These situation can be used to make estimated of the very local emissions, e.g. Hensen et al, 2000.
A.3
Inverse modeling using tall tower observations
A.3.1 Description of the COMET (inverse) transport model The COMET (CO2 MEthane Transport) model is a Lagrangian model that can be used for both predictive and inverse modeling purposes. COMET uses as the basic input backward trajectories. In the COMET model two vertical layers are distinguished, a mixing layer and a reservoir layer. The initial concentration at the start of each trajectory is taken from the averages of the calculated concentrations of a global transport model like the TM2 (Heimann, 1995) GCM. A more detailed description of the COMET model can be found in Vermeulen et al. (1997,1999).
ECN-C--04-027
39
In forward mode the COMET model retrieves per time step the emissions for the grid cells under the current circular source area and calculates the concentration changes and isotopic composition for the modeled components in the column of air with the current mixing layer height. Emission data is retrieved from emission inventories per source category on a regular grid, in order to model CO2 we use for fossil fuels emissions de EDGAR 3.2 database (Olivier et al, 2001), for sea and ocean fluxes the Takahashi dataset (Takahasi et al, 1997). Net ecosystem exchange is modelled with a simple PEM model with a superimposed daily cycle. For methane emissions the TNO METDAT database (Berdowski et al, 2000) is used. This emission database has a resolution of 0.1x0.1 degree.
Figure A4-8
The concept of the COMET model in forward mode.
The COMET Model domain is Western Europe, between –10 and 40° Eastern longitude and between 35 and 70° Northern latitude. The COMET model currently uses hourly 3D 96 hour backward trajectories (calculated with the Flextra model (Stohl et al., 1998)), the trajectories are calculated from 3D analysed meteorological fields from the ECMWF weather prediction model at 3 hour time step, a 1° x1° degree horizontal resolution and 31 vertical levels. Data is retrieved from the ECMWF MARS database on a monthly basis, using routines specifically tailored for the input files in GRIB format that the FLEXTRA model require. From the ECMWF 3D meteorological fields of windspeed, temperature and humidity a separate routine analyses the mixing layer height for each 1x1 degree gridcel using a critical Richardson number (<0.25) scheme . Then the FLEXTRA trajectory data is read and the mixing layer data and additional meteorological information is added to each trajectory record. The COMET model then uses offline from the FLEXTRA model the adapted trajectory data for each receptor point with the included meteorological information like mixing layer height, temperature, precipitation, cloud cover, global radiation, wind speed and direction at the surface. The concentration at the receptor point can be calculated directly in the models forward mode. In backward mode the models stores the contribution of aggregated source areas, the so called influence factors per trajectory. These source-response factors per trajectory and per source area form a so called source-response matrix that can be inverted in order to use it to derive best estimates of the emissions of the source areas that fit with the observations. This inversion is carried out by the SVD method. 40
ECN-C--04-027
A.3.2 Description of the SVD method Measurements are described as a function of the parameters that need to be determined by the COMET model. In our calculations no a priori information on the unknown emissions is given, and the set of equations is directly solved by (Truncated) Singular Value Decomposition (Press et al., 1987). The transport model is used to determine the contribution of each of the source areas under evaluation to the measured concentration. The measured excess of methane (X) above the background concentration (Cback) for each observed concentration (C, 1..t) is the sum of the product of all emission strengths (E) times the contributions (c) of the respective source areas (1..i), as expressed in eq. 1. This simple linear model is only valid for non-reactive and nondepositing components, otherwise unresolvable a-linearities will occur. The C matrix is delivered by the models, the X vector is derived from the measurements and the SVD technique is applied to obtain the best fit for the emission vector E between observation and model. The problem is often overdetermined, this means that conflicting information may be in the dataset and that the solution is uncertain. It could also be that emissions of certain areas contribute little to the total excess concentrations (much less than the measurement precision), but that a variation with very strong factors of the emission strengths of these areas leads to a marginally better solution. (T)SVD allows to identify those areas by objective criteria and to exclude their influence from the solution, without to much cost for the final fit between data and observations. n
X1 = C1 − Cback1t = ∑ Ei . ci ,1 i =1
M n
X t = Ct − Cback1t = ∑ Ei . ci,t
(1)
i =1
so X1 E1 c1,1 L c1,i M = M . M O M X t Ei ct ,1 L ct ,i To estimate the error in the calculations of the inverse model a statistical Monte Carlo routine is applied. In this method the model data is distorted with a user defined random factor (here 50% is assumed), corresponding with an estimated accumulated model error. This accumulated error represents the combined effects of all possible error sources. On this distorted matrix the SVD matrix inversion is applied again. This is repeated several times with newly generated random errors on the matrix elements. The standard deviation in the calculated emissions per area gives a measure of the uncertainty of the result for that area. In the inverse model the emissions are calculated for aggregate areas on a more coarse 1/2° by 1° grid (roughly equivalent to a 60 by 60 km grid). Close to the Netherlands a fine resolution is taken starting with aggregates of single cells and farther away from the Netherlands the aggregates contain more cells. The aggregate source areas also resemble in shape the expected emission or country pattern.
A.3.3 A demonstration of the method: COMET and CH4 The development of the COMET model is still ongoing. Inversions have been carried out however in the past with an earlier version of COMET and somewhat less accurate meteorological driving information, using both the methane concentration data for Cabauw in the period from 1993-1997 and concentrations in 1996 at four stations in Europe (Cabauw, Petten, Heidelberg and London). The four European stations were carefully intercalibrated (Nisbet et al, 1997). The source receptor matrix was calculated using the COMET model v1.0, using 2D trajectory data and mixing layer heights provided by EMEP and calculated with the Norwegian Weather prediction model. This data can be considered of a lower quality than the currently used ECMWF data, but was the best available data for us at that time.
ECN-C--04-027
41
The resulting emissions per country for both inversion runs are shown in Table A1. This table also shows existing emission inventories. Also the relative uncertainty in the country totals as given by the Monte Carlo error analysis is shown. When only one year of data for one station is used the solution is very unstable (not shown here). The long record from 1993-1997 at Cabauw allows for a quite stable solution, while the use of concentration data for four stations allows a fairly stable result with only one year of data. The absolute deviations between inventory data and the calculated emissions fall within the uncertainty range of the inventories, which estimated to be 30-50%. The emissions derived for the North Sea continental shelf (CS) could not be resolved using the 1 year data of the 4 stations. The transport model indicates that only 25% of the excess concentrations in Cabauw originates from the Netherlands. So the Cabauw concentration signal contains quite some information of the sources outside the Netherlands, covering an area of most of NW-Europe. Table 4
Emission estimates from inventories and inverse COMET in Tg CH4.yr-1
Country
Lotos
EDGAR
METDAT
NL D B GB EI CS
0.85 3.9 0.65 2.5 0.57 0.26
1.0 6.2 0.59 3.9 0.65 incl.
0.88 5.5 0.66 4.6 0.85 0.36
1
1
2
3
COMET Cabauw 93-97 Err% 0.76 1 3.8 1 0.67 1 5.9 1 1.5 1 0.54 3
Builtjes et al., 1992 2 Olivier et al., 1994 3Berdowski et al., 1998.
COMET 4 stations 1996 Err% 1.5 4 3.2 4 0.41 4 7.7 2 0.43 14 <0 23
The inverse model was used to derive simulated concentration data based on the relations between constructed emission data and the relations as derived by the inverse model between emission strength per emission area and the response in excess concentration at the arrival point. When this constructed concentration data is fed to the inversion routine this should lead to an exact representation of the emissions used to construct this concentration data, and this was indeed found to be the case. The Monte Carlo error analysis then allows to estimate the relation between the uncertainty in the calculated emissions and the uncertainties in concentration measurements and the model errors. As a theoretical exercise, emission fields were constructed also for nitrous oxide and HFC’s, where the nitrous oxides were assumed constant above land and freshwater at a value of the average current estimate for Europe, and zero above sea areas. The HFC emission was taken as zero over the whole area except for five point sources of different sizes in west and south-east France, Ireland, south Germany, Italy and Denmark. The five sources represent hypothetical production facilities with relatively large sources and the diffuse emission of all the other sources is thought to be relatively small per unit area and not leading to detectable changes in concentration above the base signal. With estimated precisions of the concentration measurements for these compounds and different model errors the Monte Carlo error analysis yielded the results depicted in table 3. Trajectory and mixing layer data was used for 1993-1997 for arrival point Cabauw. For methane the METDAT database was used as the a priori emission database. The error in methane concentration measurements was set at 5 ppb. The error in nitrous oxide concentrations was set at 1 ppb. The error in HFC concentration was set at 0.5%. Table 5:
Estimated maximum precision (%) of non-CO2 greenhouse gas emissions estimates using COMET inverse mode
Model error CH4 N2O HFC’s
50% 4 year 10 30 10
10% 1 year 20 unstable 30
4 year 5 15 <5
1 year 10 unstable 10
From Table 5 we can conclude that the use of longer time series allows to reduce the error in the emission estimates caused by model and measurement imperfections. Reducing the model im-
42
ECN-C--04-027
perfection from 50% to 10% reduces the uncertainties in the emissions with a factor of two. Going from 1 to 4 years of data also reduces the uncertainty with a factor of 2. The combined effect allows for theoretical uncertainties in the emission rates of about 5-10%. Uncertainties are largest for N2O because of its high background concentration compared to the size of the variation in concentration signal caused by the sources.
References Berdowski JJM, Draaijers GPJ, Janssen LHJM, Hollander JCTh, Loon M van, Roemer MGM, Vermeulen AT, Vosbeek M, Visser H., 2001. Sources, regional scaling and validation of methane emissions from the Netherlands and NorthWest Europe. Bilthoven, RIVM, report 410200084. Braswell BH, Schimel DS, Linder E & Moore BIII, 1997. The response of global terrestrial ecosystems to interannual climate variability. Science, 278, 870-872. Fan S-M., Gloor M., Mahlman J., Pacala S., Sarmiento J.L., Takahashi T., Tans P.P., A Large Terrestial Carbon Sink in North America Implied by Atmospheric and Oceanic CO2 Data and Models, Science, 282, 442-446 (1998). Feely RA, Wanninkhof R, Takahashi T & Tans P, 1999. Influence of El Nino on the equatorial Pacific contribution for atmospheric CO2 accumulation. Nature, 398, 597-601. Gloor M, Fan S, Pacala S & Sarmiento J, 1999. A model-based evaluation of inversions of atmospheric transport, using annual mean mixing ratios, as a tool to monitor fluxes of nonreactive trace substances like CO2 on a continental scale. J. Geophys. Res., 104, 14245-14260. Gloor M, Bakwin P, Hurst D, Lock L, Draxler R & Tans P, 2001. What is the concentration footprint of a tall tower? J. Geophys. Res., 106, 17831-17840. Heimann M, The global atmospheric tracer model TM2, Hamburg, Max Planck Institut für Metorologie (1995). Hensen A., Dieguez Villar A., Vermeulen A.T., Emission Estimates Based on Ambient Concentrations Measured at a 200m High Tower in the Netherlands 1995-1997. In: Ham J. van, et al (eds), Proc. 2nd Symp. on Non-CO2 Greenhouse Gases, Dordrecht, Kluwer Ac. Publ. (1999). Olivier, J.G.J. and J.J.M. Berdowski (2001) Global emissions sources and sinks. In: Berdowski, J., Guicherit, R. and B.J. Heij (eds.) ‘The Climate System’, pp. 33-78. A.A. Balkema Publishers/Swets & Zeitlinger Publishers, Lisse, The Netherlands. ISBN 90 5809 255 0. Nisbet E.H., Levin I., Wyers G.P. and Roddy A., 1997.’Quantification of the West European Methane Emissions Budget by Atmospheric Measurements’. Final report EU DGXII, EV5V-CT94-413. Press WH, Flannery BP, Teukolsky SA & Wetterling WT, 1989. Numerical Recipes in Pascal. Cambridge University Press. ISBN 0521387671. Rayner PJ, Enting IG, Francey RJ & Langenfels R, 1999. Reconstructing the recent carbon cycle from atmospheric CO2, d13C and O2/N2 observations. Tellus, 51B, 213-232. Takahashi, T., R. A. Feely, R. F. Weiss, R. H. Wanninkhof, D. W. Chipman, S. C. Sutherland, and T. T. Takahashi (1997) Global air-sea flux of CO2: An estimate based on measurements of sea-air pCO2 difference, Proceedings of the National Academy of Sciences, USA., 94, 8292-8299 Stohl A., 1998. Computation, accuracy and apllications of trajectories - a review and bibliography, Atm Env 32, 947-966. Vermeulen A.T., et al., 1997. Validation of methane emission inventories for NW-Europe. Petten, ECN, report nr ECN-C-96-088. Vermeulen A.T., Eisma R., Hensen A.,Slanina J., 1999, Transport model calculations of NW-European methane emissions. Envir. Science & Policy 2, 315-324. Vermeulen A.T., Eisma R., Hensen A., Slanina J., Transport Model Calculations of NW-European Methane Emissions. Env Sci & Policy 2, 315-324 (1999). Wofsy SC, Goulden ML, Munger JW, Fan S-M, Bakwin PS, Daube BS Jr., Bassow SL & Bazzaz FA, 1993. Net Exchange of CO2 in Midlatitude forests. Science, 260, 1314-1317.
ECN-C--04-027
43
BIJLAGE B
BIJDRAGE AAN JRC-ISPRA WORKSHOP OCT 2003
VERIFICATION OF EUROPEAN AND NATIONAL GHG INVENTORIES BY INVERSE MODELLING: STATE OF THE ART AND FURTHER RESEARCH NEEDS Peter Bergamaschi, Alex Vermeulen, Anke Herold, Frank Dentener, Klaus Radunsky and Alistair Manning
State of the art / existing studies Only during the last few years attempts have been made to derive top-down emission estimates for individual European countries by inverse modelling. This is mainly due to the big challenges to atmospheric models, in particular the requirement to simultaneously provide high spatial resolution (for resolving national boundaries) and to provide consistent global 3D tracer fields (as typically for all long-lived GHG local emissions have global influence). At the same time, continuous GHG observations over the European continent have been quite sparse, as the focus of most international GHG measurement programmes (e.g. US-American NOAA/CMDL) has been mainly to collect clean background air. This article is summarizing the material presented and discussed at the workshop, including also material available in the literature. We restrict ourselves to studies on European countries including estimates provided for the whole European continent. For an overview about global inverse modeling studies the reader is referred e.g. to the summary given by Bergamaschi and Houweling [this report, chapter 4] and the references given therein. Currently available inverse modelling studies for European countries (and the European continent) are compiled in Annex 2 and are discussed in the following separately for the 6 Kyoto gases. CO2 Inverse studies of CO2 usually prescribe fossil emissions (as well as the annually balanced part of the seasonal variation of the biosphere and ocean fluxes) in order to derive net CO2 fluxes from the biosphere or ocean. The reason for this is that fossil emissions are quite well known (assumed uncertainty a few percent), but net fluxes from the biosphere and ocean are quite uncertain. Most inverse studies listed in Annex 2 indicate that Europe’s biosphere has acted as a small net sink for CO2 during the last two decades. For exact comparison of the different studies, it is important to note that the exact definition of the European continent may be different: Most studies consider the European continent as the land mass between the Atlantic ocean and the Ural, excluding Turkey and the Mediterranean islands with a total area of 1.04 1013 m2 [Janssens et al., 2003], but other studies define Europe more broadly (e.g. including the Mediterranean islands and part of Middle East, with a total area of 1.55 1013 m2). Furthermore, the time periods are different for the different studies. The studies of [Rödenbeck et al., 2003] and [Rödenbeck, this report, chapter 6] indicate that the interannual variability of European CO2 fluxes is considerable (with values varying between -0.5 and +0.9 Pg C / yr). The discussed top-down estimates may be compared with bottom-up estimates: E.g. [Janssens et al., 2003] estimate a small and very uncertain integrated net C uptake of European forests, croplands, grasslands and peatlands of -111 ± 280 Tg C / yr. We could not identify inverse modelling studies for individual European countries. CH4 Manning [this report, chapter 6] presented estimates of national CH4 emissions from Ireland, UK, France, Germany and BENELUX. These estimates are based on a Lagrangian dispersion model (NAME model [Ryall et al., 2001]), combined with an iterative simulated annealing technique [Manning et al., 2003], and atmospheric observations at Mace Head, Ireland. In contrast to most other inverse modelling studies, this technique does not use any a priori infor-
44
ECN-C--04-027
mation about the emissions and can therefore be considered as a completely independent verification. One implicit assumption of this technique, however, is that emissions do not exhibit any seasonal variation and are constant over a period of 3-6 years. Thus, although the atmospheric signal at Mace Head is only occasionally directly influenced by the adjacent continental European countries (see also [Bergamaschi, this report, chapter 6]), this station provides significant constraints if long averaging periods are applied. The derived CH4 emissions are consistent with UNFCCC values (see Figure1, Annex 2 and Figure 4 of [Manning, this report, chapter 6]) within the given uncertainties. We note that for the discussed countries this direct comparison of top-down estimates and UNFCCC values appears justified, as the natural CH4 emissions (and sinks) are of minor importance (see e.g. [Bergamaschi, this report, chapter 6]). The inferred uncertainty of CH4 emission is ranging from ~25 % (Ireland) to ~100% (Germany), i.e. is increasing with increasing distance from Mace Head. The absolute values are slightly different compared to the previous studies [Manning et al., 2003; Ryall et al., 2001]. [Vermeulen et al., 1999] presented CH4 emission estimates for the Netherlands, Germany, UK, Belgium and France based a two-layer Lagrangian atmospheric transport model (COMET (CO2 MEthane Transport)) that describes the uptake of methane emissions in the mixing layer along a backward trajectory path and the exchange with a reservoir layer during mixing layer height changes. Using observations from Cabauw CH4 emissions for The Netherlands of 0.76 Tg CH4 / yr are derived for the 1993-1996 period. Adding three further ground based stations (London, Petten, Heidelberg) is increasing this value to 1.5 Tg CH4 / yr (for year 1996 only). For comparison, the UNFCCC numbers for Netherlands are 1.23 - 1.17 Tg CH4 / yr for the 1993-1996 period. [Roemer et al., 2000] has estimated emissions of Germany, France, Netherlands, Belgium, UK, and Ireland, based on the LOTOS model, which is a Eulerian model over the European domain (i.e. a limited domain model with boundary fields taken offline from a global model). They use observations of 4 Dutch stations and Mace Head. A striking feature in their analysis are the very high emissions from Belgium, which the authors attribute to the lack of measurement data and the possibility of errors in the used boundary conditions. The derived emissions for the other 5 countries are within a factor of 2 within the UNFCCC values. Bergamaschi [this report, chapter 6] presents an inverse study, based on an Eulerian atmospheric zoom model with 1ox1o resolution over Europe. Using a set of 10 European stations with continuous measurements (and additional global stations), the emissions of Germany, France, BENELUX, UK and Ireland are significantly constrained. The preliminary results are consistent within about 40 % with the UNFCCC numbers for these countries (second half of year 2001).
ECN-C--04-027
45
Figure 1: National estimates of CH4 emissions by inverse modelling from various studies (see legend). For comparison, the values reported to UNFCCC are also shown (solid line; assumed uncertainty 20 % [Olivier, this report, Chapter 3] indicated by gray shaded area). In addition a small correction due to natural CH4 sources is also shown (dash-dotted line), based on the inventories used in [Bergamaschi, this report, Chapter 6]. Emission from Netherlands were estimated from the reported emissions for BENELUX by applying a factor of 0.66 for those studies which reported BENELUX as total only.
N2O Similar as for CH4, [Ryall et al., 2001], [Manning et al., 2003] and [Manning, this report, chapter 6] presented N2O estimates (Figure 2, see also Annex 2). The values are close to UNFCCC numbers for the reported countries (Ireland, UK, France, Germany, BENELUX, Denmark) within the uncertainties (see error bars in Figure 2). However, N2O has also some natural sources (e.g. temperate forest and grassland soils) which need to be better quantified for a quantitative assessment.
46
ECN-C--04-027
Figure 2: National estimates of N2O emissions by inverse modelling. For comparison, the values reported to UNFCCC are also shown (solid line; assumed uncertainty 30 % [Olivier, this report, Chapter 3]: gray shaded area).
HFCs Manning [this report, chapter 6] reports national emission estimates for HFC-125, HFC134a , and HFC-152a (Figure 3, Annex 2). For HFC-125, and HFC134a there is a reasonable agreement with UNFCCC numbers for several countries within the uncertainties of top-down and bottom-up estimates, respectively. For HFC-152a, however, some discrepancies were identified: E.g. for France the top-down estimates yields HFC-152a emissions of 51-190 Tg / yr (however within uncertainties of ~100%) for 1995-2002, while only 0.7-2.5 Tg / yr are reported to UNFCCC for 1996-2001 [Gugele et al., 2003]. In contrast, top-down estimates for Germany are about a factor of ~4-5 lower than UNFCCC values. On the European scale top-down estimates for HFC-125, HFC-152a, HFC134a are also provided by [Reimann et al., 2004]. Within the given uncertainties these estimates are consistent to those of [Manning, this report, chapter 6].
ECN-C--04-027
47
Figure 3: National estimates of HFC-134a emissions by inverse modelling. For comparison, the values reported to UNFCCC are also shown (solid line; assumed uncertainty 30 %: grey shaded area).
PFCs and SF6 No inverse modelling studies have been identified for PFCs, and SF6. It is worthwhile mentioning though, that SF6 is widely used as tracer to check the transport characteristics of atmospheric models (assuming that at least the totals for the Northern and Southern hemisphere are well known to check interhemispheric exchange in atmospheric models [Denning et al., 1997]). B.1
Bottom-up inventories for use in inverse models.
In general, most inverse techniques require a priori information on the emissions totals and their spatial and temporal distribution. This is due to the fact, that typically the inverse problem is underdetermined, i.e. the numbers of unknowns (emissions per grid cell) is usually larger than the number of observations. A notable exception are the studies of [Manning, this report, chapter 6], Vermeulen [1999] and [Manning et al., 2003], where large averaging time periods are used (3-6 years) and the assumption, that over this period emissions per grid cell are constant. Observations, on the other hand, are taken on a high-frequency (e.g. 6-hourly averages), thus that finally the number of observations is larger than the number of model grid cells. An obvious advantage of this approach is, that it is indeed independent from any bottom-up inventory. Nevertheless, the mentioned assumptions may not always be valid, as eg natural CH4 sources (eg wetlands) may exhibit considerable seasonal variations. In addition, trends and interannual variations may be significant. Most other studies, however, use comparable time scales for parameters and observations (e.g. use monthly mean values in order to constrain monthly emissions) and therefore require a priori information on the sources and sinks and their uncertainties. Presently, one of the most widely 48
ECN-C--04-027
used inventory for anthropogenic sources in atmospheric (forward and inverse) models is the EDGAR database, which provides global emissions on a 1ox1 o grid (EDGAR V3.2 for years 1990 and 1995). At present, this database, does not provide any seasonal variation, which however might be important at least for some sources (e.g. CH4 from rice paddies and biomass burning). The METDAT database [Berdowski et al, 2001] contains emissions for CH4 on a 10x10 minute resolution and time factors per source category for monthly, weekday and hourly variations. UNFCCC inventories are reported only as national totals and can therefore not readily be used in atmospheric models. A first approximation, however, could be to use existing gridded inventories in order to disaggregate the UNFCCC emissions. Another limitation of the UNFCCC inventories is that they do not cover all countries Bottom-up inventories - Future requirements: The inverse modelling approaches which depend on a priori information on emissions could be improved by further improvements of the bottom-up inventories. For most inverse (and forward) models, inventories on regional scale should fulfill the following requirements: • spatial disaggregation of 1ox1o or better • temporal disaggregation for sources with significant seasonal cycles (monthly or higher temporal resolution) • in case of strong diurnal cycles this information should also be provided • uncertainty estimates, ideally also spatially (and temporally) disaggregated • Correlation of uncertainties: Most inverse approaches require the full covariance matrix (which contains the uncertainties of all model parameters and the correlations between all parameters in space and time). So far, usually very simple approximations are made, eg. assuming either no correlation between different grid cells and time steps, or a exponentially decreasing correlation with a certain spatial and / or temporal correlation length. Providing realistic covariance matrices, which consider all uncertainties and their correlations, however, will be a highly demanding task. If anthropogenic emission should be verified by inverse modelling it will be important to address also natural sources and sinks, which contribute significantly to the total budgets of CO2, CH4, and N2O. Most of these natural sources (e.g. CO2 exchange with biosphere, CH4 from wetlands, N2O from soils) are highly sensitive to climate parameters and require sophisticated biospheric models, eg [Heimann et al., 1998; Walter and Heimann, 2000]. The further developments of such models should have high priority.
B.2
Atmospheric observations
European atmospheric observations are compiled in Annex 3. For continental measurements close to sources it is very useful to measure quasi-continuously (eg. temporal resolution 1 hour or better) in order to resolve the usually strong diurnal and synoptic variations. In particular when relating the atmospheric signal to emissions, this short-term variability contains valuable information (eg. [Bergamaschi, this report, Chapter 6]) While the early atmospheric measurement programs have focused mainly on clean background air, continental measurements have been quite sparse in the past. A notable exception is e.g. the UBA / GAW station Schauinsland (Black Forest, Germany), with 30 years of continuous CO2 measurements (the longest available continental CO2 record [Schmidt et al., 2003]). Only during the last few years, efforts have been intensified to setup more European observations, in particular for CO2 under the framework of CARBOEUROPE / AEROCARB/CHIOTTO. More than 20 stations with continuous CO2 measurement are currently in operation within or close to the EU-
ECN-C--04-027
49
15 countries (see Annex 3). Most, but not all of these CO2 stations are linked to the AEROCARB project. Less stations only are available for continuous CH4 measurements (around 10 within or close to the EU-15 countries). The CH4 measurements, however, are less coordinated and many activities are part of national programs (e.g. 5 CH4 German stations are operated by the German Umweltbundesamt). Only at a few single stations measurements of N2O, HFCs, PFCs and SF6 are performed. A very promising approach is the setup of 9 tall towers throughout Europe within the EU project CHIOTTO [A. Vermeulen, this issue, Chapter 5] for continuous measurements of vertical gradients of CO2, CH4, , CO, N2O, SF6 and grab samples of relevant connected (isotopic) tracers like 13CO2, 13CH4, O2/N2.Furthermore, the monitoring of GHGs is receiving more attention, e.g in the framework of GMES and the further development of the Global Climate Observing System. Nevertheless, the organization of GHG monitoring in Europe is still very heterogeneous: Within Germany, the German Umweltbundesamt is maintaining an operational network for GHG measurement. Most other monitoring activities in Europe, however, are part of European and national research projects. In additional, some stations are part of international networks as GAW and AGAGE. This heterogeneous structure mainly implies three problems: - Most European research projects currently have very restrictive data policies. Therefore, data collected within these projects are often not available to scientists outside these projects
Atmospheric observations - Future requirements: •
Ensure further operation of quasi-continuous (e.g. hourly) measurements of CO2, CH4, N2O, HFCs , PFCs and SF6 mixing ratios at existing stations. For use in inverse models (in particular for trend assessments), long-term data sets (avoiding any data gaps) are very important. Changing combinations of atmospheric monitoring sites will introduce biases in the analyses and should be avoided.
•
As many regions in Europe are still poorly covered (in particular South / East Europe, but also Northern Europe) further measurement stations need to be setup. This is in particular true for all non-CO2 GHGs. For an European operational verification system of national emission inventories a minimum of ~1-5 stations / country (depending on country size) would be required. Also, as every station has its own characteristics (e.g. topography, local wind systems which may be difficult to model) a certain redundancy of stations would be desirable.
•
In order to ensure high-quality long-term measurements (including careful calibration and intercomparison) on the medium-term a European network/framework for monitoring of GHGs should be setup. Such a network could e.g. be comparable to the EMEP network for verification of the Convention on Long Range Transboundary Air Pollution (LRTAP), and should integrate the monitoring activities of the current GHG research projects and networks (as UBA, GAW, AGAGE, EMEP, CARBOEUROPE IP). Such a European framework could strongly benefit from synergies, using existing infrastructure in an optimal and most cost-effective way.
•
Further research on remote sensing has high priority. In particular for CH4 and CO2 remote sensing may provide significant constraints if observation with accuracies of better than 1 % can be provided. Beside the accuracy requirement, however, also the global coverage is very important. At the moment, clouds still represent a major problem, limiting the number of useful pixels typically to a small percentage of the whole data set
50
ECN-C--04-027
- Although considerable efforts have been undertaken in intercomparisons between different stations / laboratories (e.g. within the CARBOEUROPE projects for CO2), full intercomparability of all European GHG measurements (e.g. use of same calibration scale) may not always be guaranteed. - The various projects and institutions decide about their monitoring sites depending on their priorities and resources. Therefore, from a European point of view, the overall set of stations may be not the optimal one (and long-term maintenance of stations not always guaranteed) Additional information about the relative contribution of different sources (and sinks) can also be provided by isotopic measurements, utilizing the fact that for many gases certain sources (or groups of sources) have distinct isotopic fingerprints. E.g. 13C/12C and D/H ratios in methane are significantly different for CH4 from different production processes: biogenic sources (e.g. wetlands, rice paddies, ruminants), thermogenic formation (natural gas), and incomplete combustion of biomass [Wahlen, 1993]. This information can be utilised for identifying the contribution of the different source categories (e.g. [Bergamaschi et al., 1998]). For CO2, 14CO2 allows the identification of fossil sources [Levin and Kromer, 1997]. For N2O, 15N/14N and 18O/16O ratios may be used to identify anthropogenic sources [Röckmann et al., 2003]. In general, however, isotopic measurements require extremely high precision and accuracy, are usually very labor-intensive and expensive, which has restricted isotopic measurement to very few global stations (and measurement campaigns) only. The development of new techniques for automated sample processing may improve the situation [Miller et al., 2002]. At this workshop, however, isotope measurements were not covered further in detail. B.3
Atmospheric models
Various model approaches have been used for the ‘regional’ inverse modelling B.4 Lagrangian back trajectory models This approach calculates the back trajectory path for a given measurement site. In the simplest case, the back trajectories are calculated on isentropic surfaces only. More sophisticated approaches provide a full 3 D picture, including parameterised subgrid scale processes (turbulent vertical mixing, convective clouds). An obvious advantage of the back trajectory approach is that it directly provides the source attribution maps for a given monitoring site, i.e. no further differentiation techniques (such as an adjoint model) are required. Examples are the NAME model [Ryall et al., 1998] used by Manning [this report, chapter 6] and the COMET model applied in the study of [Vermeulen et al., 1999]. These models have horizontal resolutions in the order of 0.5 o -1o, i.e. resolve the national boundaries of larger countries (e.g. France) reasonably well. Typical integration times are several hours to several days, i.e. only the very recent history of an air parcel is calculated. The technique requires the determination of baseline concentrations (which are subtracted from the total signal). B.5 Lagrangian Particle Dispersion Models Particle modeling is the Lagrangian approach towards the simulation of atmospheric dispersion. Lagrangian particle models compute trajectories of a large number of individual so-called particles (not necessarily representing real particles, but small air parcels) to describe the transport and diffusion of substances in the atmosphere. Usually they are used to simulate the dispersion of air pollutants, released for instance due to an accident in a nuclear power plant. A main advantage of Lagrangian models is that there is no artificial numerical diffusion as occurring in Eulerian models. This is of special importance in the vicinity of the source of an air pollutant, where in Eulerian models the pollutant is instantaneously mixed over at least one grid box, which can cause large subsequent transport errors. The description of turbulence with parECN-C--04-027
51
ticle models is of similar quality as large eddy simulations, and it is more accurate than with Eulerian models. Particle models are independent of a computational grid and, in principle, provide infinitesimally small resolution. Practically, however, their effective resolution is determined by the number of particles released and by the resolution of meteorological input data. Tens or hundreds of thousand particles may be required to accurately describe a multiday transport episode. This posts high demands on computer resources and explains why particle models have only become popular in recent years since powerful computers became available. The basis for most of nowadays particle models was developed in a paper by Thomson (1987). An example LPDM is the FLEXPART model (Stohl et al, 1998). FLEXPART may be used in forward mode to simulate the dispersion of pollutants released at a location, or in backward mode to determine the ``influence area’’ of a given receptor point. Forward trajectories are calculated in the forward mode, backward trajectories in the backward mode. Deposition processes can be simulated in both forward and backward mode. In its latest version FLEXPART has an extensive algorithm to deal with the influence of convection in the vertical particle distribution. LDPM are much more demanding than Lagrangian trajectory models and most Eulerian models. Main of advantage of LDPMs is that the output of the models can be converted to any resolution of output grid in horizontal or vertical direction, better matching the usually pointwise receptor points of concentration measurements, as long as the statistics of the number of particles per grid cell allow for significant certainty, which can be easily extended by modelling more particles (to the expense of computing overhead, typically hundreds of thousands to millions of particles are used in a simulation). Because of the linear nature of the models, the output of LDPMs like FLEXTRA can be relatively easy used to solve inverse model exercises. Another advantage is take the nature of the descriptions of the transport allows to track for each particle from what source (type/category/place it originated) The added statistics of the number of particles per output grid cell allows for statistical evaluation of the model errors and uncertainty of the predicted concentrations. The following approaches are based on Eulerian models1 : B.6 Coupled mesoscale / global models Due to the computational requirements models of very high spatial resolution are usually applied on a limited area only. In particular for simulations of longlived GHGs, however, the global distribution plays a significant role. Therefore, the application of mesoscale / limited area models requires the provision of global 3D tracer fields from global models. Such an offline coupling, however, ignores any feedback from the mesoscale model domain to the global domain. Also, when applied in inverse mode, such a coupling is very challenging in terms of full consistency between mesoscale and global domain. Examples for the application of this approach are the studies of [Chevillard et al., 2001; Vermeulen et al., 1999]. B.7 Zoom models In order to overcome the problems of offline coupling zoom models have been developed. In this approach the region of interest is simulated at high spatial resolution, embedded into a global model, run at a lower resolution. An example is the TM5 model [Berkvens et al., 1999; Krol et al., 2004] (see also TM5 web site http://www.phys.uu.nl/~tm5/). This approach guarantees consistent boundary conditions for the zoom region and the two-way nesting accounts for feedback effects from the zoom region to the global domain (which is in particular important at the zoom region boundaries, with frequent in- and outflow of airmasses). Another recent zoom model is the LMDZ model (http://www.lmd.jussieu.fr/~lmdz/). In contrast
1
Lagrangian approach: Air parcel (or particle) moves with local wind; in order to treat dispersion and sub grid scale mixing processes larger number of individual particles may be used. Eulerian approach: model simulates species concentration in a fixed 3D model grid
52
ECN-C--04-027
to the TM5 models, the LMDZ models has only one model domain, with a dynamically changing resolution. B.8 High-resolution global models With continuously fast development of computer techniques it will also be an option to run high resolution models (e.g. 1o x 1o or better) on a global domain. It has to be kept in mind, however, that inverse calcuations typically require 1-2 orders of magnitude more CPU resources compared to forward simulations. Nevertheless this global high resolution approach is feasible, with the advantage of providing high resolution global results and avoiding the zooming algorithms (which introduce an additional complexity in atmospheric models). Atmospheric models - Future requirements: •
although a spatial resolution of 1o x 1o already reasonably resolves national boundaries of larger countries, further increase of spatial resolution is desirable, in particular to better resolve local meteorological effects. ECMWF analyses (used as input for most offline transport models) are currently provided at a resolution of ~0.5o x 0.5o, but will soon be provided at higher resolution of ~0.2o x 0.2o. The resolution of the meteorological input fields limits the maximum potential resolution of offline transport models. An alternative of using offline models could be the application of full GCMs (Global Circulation Models), however implying much higher CPU requirements and the need of coupling to global meteorological analyses (such as ECMWF).
•
Recent model intercomparisons, based on 222Rn or SF6 (eg. within the EVERGREEN project) indicate that many atmospheric models simulate advective transport reasonably well. However, different models show significant differences in vertical transport. In particular further improvement of subgrid scale processes should have high priority.
•
A notorious problem in inverse modelling is the exact specification of model errors / uncertainties. Therefore, ensemble models runs (e.g. using the same model with different parameterisations, or different models for the same problem) could provide better model uncertainty estimates [Gurney et al., 2002]
Inversion techniques - Future requirements: •
The development of sophisticated data assimilation techniques (e.g. similar to those used in numerical weather prediction) should have high priority.
•
So far most inverse studies use rather crude approximations of parameter uncertainties (covariance matrix). However, these uncertainties have significant influence on a posteriori results (except for those studies, which do not use any a priori information). A particular weak point are potential temporal and spatial correlations (non-diagonal elements of covariance matrix).
B.9
Inversion techniques
The principal inversion techniques applied in regional studies are the same or similar to those in global studies (see eg. overview by [Bergamaschi and Houweling, this issues, Chapter 4]). For regional studies usually a higher temporal resolution should be used (in order to fully exploit the information content of synoptic variations). At the same time the number of parameters to be derived (emission per grid cells and time) is increasing with increasing spatial resolution. Thus the increasing numbers of both parameters and observations represents a major challenge for future inversion techniques. ECN-C--04-027
53
References Bergamaschi, P., C.A.M. Brenninkmeijer, M. Hahn, T. Röckmann, D.H. Scharffe, P.J. Crutzen, N.F. Elansky, I.B. Belikov, N.B.A. Trivett, and D.E.J. Worthy, Isotope analysis based source identification for atmospheric CH4 and CO sampled across Russia using the Trans-Siberian railroad, J. Geophys. Res., 103 (D7), 8227-8235, 1998. Berkvens, P.J.F., M.A. Botchev, W.M. Lioen, and J.G. Verwer, A zooming technique for wind transport of air pollution, pp. 31, Centrum voor Wiskunde en Informatica, 1999. Chevillard, A., U. Karstens, P. Ciais, S. Lafont, and M. Heimann, Simulation of atmospheric CO2 over Europe and Siberia using the regional scale model REMO, Tellus, submitted, 2001. Denning, A.S., M. Holzer, K.R. Gurney, M. Heimann, R.M. Law, P.J. Rayner, I.Y. Fung, S.-M. Fan, S. Taguchi, P. Friedlingstein, Y. Balkanski, M. Maiss, and I. Levin, Threedimensional transport and concentration of SF6: A model intercomparison study (TransCom 2), Tellus, submitted, 1997. Gugele, B., K. Huttunen, and M. Ritter, Annual European Community greenhouse gas inventory 1990-2001 and inventory report 2003, European Environment Agency, Copenhagen, 2003. Gurney, K.R., R.M. Law, A.S. Denning, P.J. Rayner, D. Baker, P. Bousquet, L. Bruhwiler, Y.H. Chen, P. Ciais, S.-M. Fan, I.Y. Fung, M. Gloor, M. Heimann, K. Higuchi, J. John, T. Makl, S. Maksyutov, K. Masarie, P. Peylin, M. Prather, B.C. Pak, J. Randerson, J. Sarmiento, S. Taguchi, T. Takahashi, and C.-W. Yuen, Towards robust regional estimates of CO2 sources and sinks using atmospheric transport models, Nature, 415, 626-630, 2002. Heimann, M., G. Esser, A. Haxeltine, J. Kaduk, D.W. Kicklighter, W. Knorr, G.H. Kohlmaier, A.D. McGuire, J. Melillo, B. Moore III, R.D. Otto, I.C. Prentice, W. Sauf, A. Schloss, S. Sitch, U. Wittenberg, and G. Würth, Evaluation of terrestrial carbon cycle models through simulations of the seasonal cycle of atmospheric CO2: First results of a model intercomparison study, Global Biogeochem. Cycles, 12 (1), 1-24, 1998. Janssens, I.A., A. Freibauer, P. Ciais, P. Smith, G.-J. Nabuurs, G. Folberth, B. Schlamadinger, R.W.A. Hutjes, R. Ceulemans, E.-D. Schulze, R. Valentini, and H. Dolman, Europe’s terrestrial biosphere absorbs 7-12 % of European anthropogenic CO2 emissions, Science Express, 300 (5623), 2003. Krol, M.C., S. Houweling, B. Bregman, M. van den Broek, A. Segers, W. Peters, F. Dentener, and P. Bergamaschi, A global chemistry transport model that allows two-way nested zooming: algorithm and applications, Atmos. Chem. Phys. Disc., to be submitted, 2004. Levin, I., and B. Kromer, Twenty years of atmospheric 14 CO2 observations at Schauinsland station, Germany, Radiocarbon, 39, 205-218, 1997. Manning, A.J., D.B. Ryall, R.G. Derwent, P.G. Simmonds, and S. O’Doherty, Estimating European emissions of ozone-depleting and greenhouse gases using observations and a modeling back-attribution technique, J. Geophys. Res., 108 (D14), 4405, doi:10.1029/2002JD002312, 2003. Miller, J.B., K.A. Mack, R. Dissly, J.W.C. White, E.J. Dlugokencky, and P.P. Tans, Development of analytical methods and measurements of 13 C/ 12 C in atmospheric CH 4 from the NOAA Climate Monitoring and Diagnostics Laboratory Global Air Sampling Network, J. Geophys. Res., 107 (D13), doi:10.1029/2001JD000630, 2002. Reimann, S., D. Schaub, K. Stemmler, D. Folini, M. Hill, P. Hofer, B. Buchmann, P.G. Simmonds, B.R. Greally, and S. O’Doherty, Halogenated Greenhouse Gases at the Swiss High Alpine Site of Jungfraujoch (3580 m asl): Continuous Measurements and their Use for Regional European Source Allocation, J. Geophys. Res., in press, 2004. Röckmann, T., J. Kaiser, and C.A.M. Brenninkmeijer, The isotopic fingerprint of the preindustrial and the anthropogenic N2O source, Atmos. Chem. Phys., 3, 315-323, 2003, 3, 315-323, 2003.
54
ECN-C--04-027
Rödenbeck, C., S. Houweling, M. Gloor, and M. Heimann, CO2 flux history 1982-2001 inferred from atmospheric data using a global inversion of atmospheric transport, Atmos. Chem. Phys., 3, 1919-1964, 2003. Roemer, M., M. van Loon, and G. Boersen, Methan emission verification on a national scale, TNO, Apeldoorn, The Netherlands, 2000. Ryall, D.B., R.G. Derwent, A.J. Manning, P.G. Simmonds, and S. O’Doherty, Estimating source regions of European emissions of trace gases from observations at Mace Head, Atmospheric Environment, 35, 2507-2523, 2001. Ryall, D.B., R.H. Maryon, R.G. Derwent, and P.G. Simmonds, Modelling long-range transport of CFCs to Mace Head, Q.J.R. Meteorol.Soc., 124, 417-446, 1998. Schmidt, M., R. Graul, H. Sartorius, and I. Levin, The Schauinsland CO2 record: 30 years of continental observations and their implications for the variability of the European CO2 budget, J. Geophys. Res., 108 (D19), 4619, doi:10.1029/2002JD003085, 2003. Vermeulen, A., R. Eisma, A. Hensen, and J. Slanina, Transport model calculations of NW European methane emissions, Env. Sci. & Policy, 2, 315-324, 1999. Wahlen, M., The global methane cycle, Annu. Rev. Earth Planet Sci., 21, 407-426, 1993. Walter, B.P., and M. Heimann, A process-based, climate-sensitive model to derive methane emissions from natural wetlands: Application to five wetland sites, sensitivity to model parameters, and climate, Global Biogeochemical Cycles, 14 (3), 745-765, 2000.
ECN-C--04-027
55
