Observation-based estimates of fossil fuel-derived CO2 emissions in the Netherlands using Δ14C, CO and 222Radon

Surface emissions of CO2 from fossil fuel combustion (FFCO2) are estimated for the Netherlands for the period of May 2006–June 2009 using ambient atmospheric observations taken at station Lutjewad in the Netherlands (6◦21E, 53◦24N, 1 m. a.s.l.). Measurements of 14C on 2-weekly integrations of CO2 and CO mixing ratios are combined to construct a quasi-continuous proxy record (ΦFFCO2 ∗ ) from which surface fluxes (ΦFFCO2 ∗ ) are determined using the 222Rn flux method. The trajectories of the air masses are analysed to determine emissions, which are representative for the Netherlands. We compared our observationally based estimates to the national inventories and we evaluated our methodology using the regional atmospheric transport model REMO. Based on 3 yr of observations we find annual mean ΦFFCO2∗ emissions of (4.7 ± 1.6) kt km−2 a−1 which is in very good agreement with the Dutch inventories of (4.5 ± 0.2) kt km−2 a−1 (average of 2006–2008).


Introduction
Human activities have induced global warming by emitting large amounts of the long-lived greenhouse gases carbon dioxide (CO 2 ), methane (CH 4 ), nitrous oxide (N 2 O) and halocarbons. From this group, CO 2 is the largest contributor to global warming (Forster et al., 2007). To prevent global warming to rise to dangerous levels, most countries have joined an international treaty, the United Nations Framework Convention on Climate Change (UNFCCC), with the intention to reduce their greenhouse gas emissions. Currently, 184 parties of this convention have also signed the Kyoto protocol subjecting themselves to legally binding targets to reduce their emissions by on average 5% by the year 2012 (calculated as the average of 2008-2012) compared to 1990. At the 15th Conference of Parties in Copenhagen in 2009 much larger reduction targets were negotiated but a legally binding agreement could not yet be agreed upon. Nevertheless, the Netherlands declared to reduce its emissions by 30% by the year 2020 (compared to 1990) and the European Union has agreed to a reduction target of 20-30%.
Currently, these reduction targets are not validated by means of an independent approach. Parties only report their annual greenhouse gas emissions based on statistical inventories (UNFCCC, 2009), by estimating the emission or uptake of various source and sink processes, and adding them all up. However, these bottom-up statistics can potentially be biased if the emission or uptake of relevant sources or sinks are incorrect or missing. Furthermore, the uncertainties can be as large as, or even larger than, the reduction target itself (Rypdal and Winiwarter, 2001;Levin et al., 2003). In principle, the only truly independent method for verification of emission reductions is by observing the changes in the atmosphere. However, observations of CO 2 mixing ratios alone are not sufficient because of the very large natural CO 2 fluxes and secondly, the atmosphere integrates all emissions in space and time. This makes it difficult to accurately determine the magnitude of the fluxes (from surface to atmosphere or vice versa) and to aggregate them over a certain area. One commonly applied method for estimating surface fluxes from ambient mixing ratios is to combine ambient measurements with atmospheric transport models; the so-called inverse modelling method (Rayner et al., 1999;Bousquet et al., 2000;Rödenbeck et al., 2003;Bergamaschi et al., 2005;Peylin et al., 2005;Baker et al., 2006). In spite of the mathematical elegance of this method there are still large uncertainties involved regarding (e.g.) estimations of the boundary layer height, modelled transport in the atmosphere and errors related to the resolution of the model (Engelen et al., 2002;Rödenbeck et al., 2006;Tolk et al., 2008). Furthermore, the method is highly sensitive to measurement biases between different observation stations as these would be translated by the model into a very strong source or sink in-between. This method is therefore not (yet) suitable and is currently not applied for estimating the fossil fuel based CO 2 .
A more suitable observation-based method for estimating surface fluxes is the 222 Radon ( 222 Rn) flux method (Levin, 1984;Thom et al., 1993;Schmidt et al., 1996). 222 Rn is a radioactive noble gas (its half-life time is 3.8 d) that is produced at a constant rate from 226 Radium, which is relatively uniformly distributed in all soils. When released into the atmosphere, 222 Rn experiences the same atmospheric circumstances (i.e. transport and dilution through mixing) as all other gases released from the surface. If the Rn flux is known and its atmospheric concentration measured, the ratio between those two can be determined, and subsequently applied to calculate surface fluxes from mixing ratios of other gases in the atmosphere. The main condition for using this method is that the 222 Rn soil emission has to be well known in time and space. This method has been applied successfully to estimate the emissions of CH 4 and N 2 O in the Netherlands (Van der Laan et al., 2009b). In this paper, we will demonstrate the method for CO 2 emissions from fossil fuel combustion.
After a description of our measurement site Lutjewad and the equipment used (Section 2), we will describe the use of the radioactive isotope 14 C to identify the CO 2 from fossil fuel combustion (FFCO 2 ) in Section 3.1. Since CO 2 from fossil fuels contains no 14 C anymore due to its high age, 14 C can be used as a proxy for fossil fuel CO 2 emissions (Levin, 1987;Zondervan and Meijer, 1996;Turnbull et al., 2006). Because 14 C is too labour intensive to allow continuous observations, we apply a second tracer: Carbon Monoxide (CO). CO has proven to be a valuable proxy for FFCO 2 Levin et al., 2003;Gamnitzer et al., 2006;Turnbull et al., 2006) as its sources are very closely linked to that of FFCO 2 . Any hydrocarbon oxidation process with CO 2 as an end product is to some extent associated with CO production (Gamnitzer et al., 2006). We therefore calibrate our carbon monoxide (CO) mixing ratios to FFCO 2 , which are integrated samples over 2 weeks, and construct a high-resolution proxy for fossil fuel derived CO 2 (FFCO 2 * ). We subsequently apply the 222 Rn flux method to calculate the surface emissions FFCO 2 * (Section 3.2) and make a distinction between Dutch emissions and emissions from abroad by looking at their trajectories. Similarly, we distinguish between emissions from the northern part of the Netherlands (where station Lutjewad is located) and emissions from the densely populated centre and southern part from where we expect the highest emissions. For comparison, and for testing the methodology, the same exercise was performed on modelled concentrations of 222 Rn and FFCO 2 from the regional transport model REMO (Section 3.3). Finally, the results are compared to the national inventories (Section 4).

Measurement site Lutjewad and applied instrumentation
All measurements were performed at our station Lutjewad (53 • 24 18 N, 6 • 21 13 E, 1 m above sea level). The station is located in the north of the Netherlands on the Dutch North Sea coast, about 30 km to the northwest of the city of Groningen. Measurements include quasi-continuous observations of CO 2 , CH 4 , N 2 O, SF 6 , CO and 222 Radon, automated flask sampling ( 13 C, 14 C, O 2 /N 2 , CO 2 , CH 4 , CO), Eddy covariance (CO 2 , H 2 O) and basic meteorological properties (air temperature, humidity, atmospheric pressure, wind-speed and direction and solar radiation). See also Neubert et al. (2004). Wind conditions at the site are such that most air is sampled from wind directions between southwest to west. During the period of May 2006-May 2009 the prevailing wind direction (31%) was between 195 • and 255 • and most wind speeds (35%) were between 6 and 9 m s −1 , based on observations at 60 m above ground level (the same as the sample intake) ( Van der Laan et al., 2009b). The majority of the sampled air in Lutjewad is therefore highly influenced by emissions from the Netherlands.
From the intake at a height of 60 m, air is flushed down continuously to a laboratory where the air is dried to a dew point of −50 • C, after which the analyses take place and the air is collected in flasks.
Following a travel time of about 4 min the sample air is fed to a modified Agilent HP 6890N Gas Chromatograph (GC) where separation and analyses of CO 2 and CO (also CH 4 , N 2 O and SF 6 ) take place. About six measurements are performed in 1 h and the obtained measurement uncertainty is about ±0.04-0.06 ppm for CO 2 and 0.8-1.8 ppb for CO (Van der Laan et al., 2009a).
Bi-weekly integrated samples of CO 2 are collected by absorption of CO 2 in CO 2 -free sodium hydroxide (NaOH) solutions. 14 C is then determined by conventional 14 C analysis, or in the case of small samples due to unfavourable wind conditions, in duplicate with Accelerator Mass Spectrometry (AMS) in our radiocarbon laboratory in Groningen (Stuiver and Polach, 1977;Mook and Van der Plicht, 1999). For this study, samples were collected only when wind direction was between 100 • and 250 • in order to sample continentally influenced emissions. The measurement precision is typically ±2-3 for both methods. Ambient concentrations of 222 Rn are measured (half hourly averages) with an ANSTO dual-flow loop two-filter detector (Whittlestone and Zahorowski, 1998). The total travel time from the inlet in the mast, at 60 m height, is about 10 min. The detector is a non-energy selective alpha particle counter and would detect ambient 220 Rn (half-life time of 55.6 s) as well as 222 Rn. However, this is prevented because of the relatively long travelling time from the tower inlet to the detector compared to the decay of 220 Rn, that is, roughly 10 half life times. A filter in front of the detector removes aerosols and (radioactive) decay products. The 222 Rn decay products are sampled on a second filter in a 1500 litre delay chamber, where their decays are counted by a photo-multiplier. The combined measurement uncertainty depends on the total decay counts and the uncertainty of the 222 Rn source with which the device is calibrated, and is typically 5%.

Constructing a quasi-continuous proxy for CO 2 from fossil fuels
Our first step is to identify the 2-weekly integrated CO 2 from fossil fuels (FFCO 2 ) by combining our semi-continuous CO 2 mixing ratios with the 2-weekly integrated 14 C measurements. We subsequently calibrate carbon monoxide (CO) to FFCO 2 to obtain a high-resolution proxy for FFCO 2 (FFCO 2 * ), which is needed for estimating the Dutch surface emissions of CO 2 from fossil fuels. Ambient observed mixing ratios of CO 2 (CO 2obs ) consist of a background component (CO 2bg ), a fossil fuel component (FFCO 2 ), a biosphere component (CO 2bio ) and an oceanic component (CO 2oc ). In our case we include CO 2oc in CO 2bg since we assume the 14 C gradient between Lutjewad and Jungfraujoch (which we use as a background reference station) is negligible. This gives CO 2obs = CO 2bg + CO 2bio + FFCO 2 . (1) Expressed in 14 C, which is the ratio of 14 C/ 12 C per mille deviations from the standard pre-industrial atmosphere and corrected for radioactive decay and mass dependent fractionation (Stuiver and Polach, 1977), gives CO 2obs ( 14 C obs + 1) = CO 2bg ( 14 C bg + 1) + CO 2bio ( 14 C bio + 1) where the bio-term includes both respiration and photosynthesis. As the 14 C in fossil fuels is already decayed, 14 C ff = -1 according to its definition, and thus the final term on the right is zero. This leads to Following Meijer et al. (1996) and Zondervan and Meijer (1996) we set CO 2bio = CO 2bg in eq. (3) which gives FFCO 2 = CO 2obs · 14 C obs − 14 C bg The biosphere component (CO 2bio ) consists of about 50% of autotrophic respiration, which is in close equilibrium with the background (Levin et al., 2008). The heterotrophic respiration should in principle be taken into account for terrestrial sites (Randerson et al., 2002) but is due to the large heterogeneity of the biosphere difficult to estimate precisely (e.g. Levin et al., 2003;Gamnitzer et al., 2006;Graven et al., 2009). Based on a mean terrestrial residence time of 10 yr Turnbull et al. (2006) estimated for the northern hemisphere that FFCO 2 would be underestimated due to this biospheric 14 CO 2 source between 0.2 ppm in winter to 0.5 ppm in summer. In this study we therefore add a harmonic regression fit between 0.2 ppm in winter and 0.5 ppm in summer to our FFCO 2 values from eq. 4 (ε resp ) as an approximation for the seasonal 14 CO 2bio effect on FFCO 2 .
Because 14 CO 2 is constantly produced in the upper atmosphere, the ideal background reference 14 C bg would be sampled at a station low enough to avoid variations caused by changes in the production rate (e.g. changes in cosmic radiation intensity) but high enough from the earth's surface to avoid depletion of its mixing ratio by fossil fuel CO 2 emissions (Hesshaimer, 1997;Randerson et al., 2002). The High Alpine site Jungfraujoch (JFJ) (3450 m a.s.l.) in the Swiss Alps (Levin et al., 2003;Levin et al., 2008) is about the best practical realisation available of this ideal situation, so we used the 14 C observations of JFJ as a background reference fitted with a linear trend fit with a 1-harmonic seasonality.
The 2-weekly integrated 14 C mixing ratios from our NaOHbased CO 2 sampler at Lutjewad, as well as the JFJ 14 C measurements for the period of May 2006-May 2009 are shown in Fig. 1a together with the Dutch monthly mean temperature (source: www.KNMI.nl). Very low 14 C values are observed in the (very cold) winter of 2008-2009, compared to the previous years. An extremely low 14 C value of −6.4% was measured in the first 2 weeks of January 2009 indicating high influence of fossil fuel burning and low atmospheric mixing. Fig. 1b shows the corresponding 2-weekly averaged FFCO 2 mixing ratios as calculated with eq. (4) and corrected for 14 C bio . The regression fit, a LOESS weighted smooth with a smoothing factor set to 1 yr, illustrates the clear seasonal pattern but does not capture the very low and high excursions. High FFCO 2 mixing ratios are generally observed in the winter and lower values in the summer, both influenced by atmospheric conditions (i.e. transport and planetary boundary height) as well as anthropogenic activity. The very low 14 C excursion of −6.4% in Fig 1a is shown as a very high FFCO 2 mixing ratio of 20.6 ppm. The seasonal amplitudes vary strongly between the years with about 4 ppm in 2006 to about 9 ppm in 2009. The lower values are attributed to the fact that the year 2006 was a very warm year; the winter of 2006-2007 was the warmest in the Netherlands since at least 300 yr (KNMI, 2009). The uncertainty in the FFCO 2 observations is typically 1.4 ppm. The corresponding relative error is around 20% for most observations, but can become much larger (i.e. >100%) when the 14 C observations are close to the JFJ reference (i.e. very low FFCO 2 ).
For practical reasons ( 14 C measurements are too costly and time-and labour-intensive to allow high resolution) our 14 C obs (and thus FFCO 2 ) are limited to 2-weekly integrated samples. Therefore, we calibrated our (2-weekly integrated) high resolution CO mixing ratios to 2-weekly integrated FFCO 2 and convert our CO data set to a proxy for FFCO 2 which is defined here as: FFCO 2 * . First, we determined a background for CO by fitting the daily minimum values with a weighted mean, and subtracted it from the CO mixing ratios in order to determine the CO enhancements from the background (dCO). Fig. 1c shows the 2-weekly integrated dCO observations, fitted with a weighted linear trend fit with a harmonic seasonality to guide the eye. The high excursions from the background value, which are observed in the winters, are not well represented by the fit. For each 2 week-integrated sample 14 C obs we calculated a ratio dCO/FFCO 2 which is shown in Fig. 1d. Most observed values are between 5 and 15 ppb dCO/FFCO 2 ppm. This low ratio is related to the high share of natural gas in fossil fuel consumption in the Netherlands  which is used for a major part of the electricity production and for almost all heat-ing of buildings. The seasonal cycle of dCO/FFCO 2 is therefore mainly temperature driven. In the winter periods, the fossil fuel emissions in the Netherlands are dominated by the heating of buildings which decreases the ratio dCO/FFCO 2 because very little carbon monoxide is formed when combusting natural gas. Because of the high variability and often large error bars, the data were fitted with a weighted linear trend fit with a harmonic seasonality which was used to convert the CO mixing ratios to FFCO 2 * . This fit is mainly determined by observations from mid 2007-mid 2009, as there are fewer observations in the first year due to technical difficulties. The fit is also in lesser extent influenced by the observations in the summer where the uncertainty is generally very large. Figure 1e shows the high resolution mixing ratios of FFCO 2 * which also have a clearly seasonal cycle with higher mixing ratios in the winter periods and lower values in the summers. The 222 Rn observations for the corresponding period are shown in Fig. 1f. Higher concentrations are generally observed in the winters when the atmosphere is more stable.

Event selection
To calculate surface fluxes from FFCO 2 * ( FFCO 2 * ) and distinguish between different source areas, we selected the data based on specific events according to the build-up of 222 Rn in the atmosphere which is common with stable atmospheric conditions. Figure 2 shows an example of some of these events, indicated with shaded areas, observed at our station from November 4 to November 24 in the year 2007. We selected an event when a significant departure of 222 Rn from the background concentrations was observed for at least four consecutive hours. Due to lack of vertical mixing with free tropospheric air, all surface emissions (e.g. 222 Rn and FFCO 2 ) are trapped in the lower atmospheric boundary layer. An event terminates when vertical mixing is reestablished. Events which start already before background level is regained (indicating a long period of continental influence, indicated with the number 2 in Fig. 2) are treated similarly to the other events after defining a new baseline (illustrated by the dashed lines). Similarly to Van der Laan et al. (2009b) we assume that the meteorological circumstances (e.g. wind speed) remain stable and the transit time is the same for the whole air mass, and use the length of an event as an indicator for the area of influence. We analysed the trajectories of the air masses with the Hysplit 4 lagrangian back trajectory model (Draxler and Rolph, 2003) using hourly Global and hourly Data Assimilation System (GDAS) metrological data (downloadable from http://ready.arl.noaa.gov/gdas1.php) for each event correspond-ing to the total duration of the event. For example, if an event sustained for 8 h we calculated a trajectory 8 h back in time. As a starting point we used the last data point of an event (before vertical mixing was re-established) as well as a point in the middle of an event in order to validate steady-state conditions during travel. Events were accepted for further analysis if the trajectories indicated that the track of the air mass was mainly (>70%) over the Netherlands in order to select events, which were dominantly influenced by Dutch emissions. Figure 3 shows a contour plot constructed from the hourly points on the trajectories of all events which were selected for further analysis. The values are normalized towards the end points of the trajectories (station Lutjewad, indicated with an "x") in order to show the relative sampling density over the footprint. Each event represents the average surface emissions along its trajectory, while the effective capture range around such a trajectory gets broader the further one goes back in time. Therefore, to determine the mean emissions of the Netherlands we would only have to collect a series of events that together represent the total area of the Netherlands, implying that such trajectories would have to start close to the borders of the country. Attributing the same weight to trajectories that start close(r) to our station as to those further reaching ones would lead to an overrepresentation of the region close to the station, which is also clear from Fig. 3. As a practical solution, we decided to divide the country into two sectors: the first relatively close (∼100 km radius) to the station (sector 1 in Fig. 3), and the second comprising the rest of the Netherlands (sector 2 in Fig. 3). We then assume the trajectories starting in the outer sector to be representative for the Netherlands as a whole, whereas those starting in the inner sector would represent the regional emissions of the north of the Netherlands. The emissions from sector 2 itself are then calculated in a straightforward way from them. The two sectors are also representative for the population density in the Netherlands with sector 1 having less than 250 inhabitants km −2 whereas the average of the whole country is about 500 inhabitants km −2 . Once this separation has been made, the calculation of the flux is just the arithmetic mean of the flux results for the individual events. Accordingly, the uncertainty is the square root of the quadratic sum of individual event uncertainties, divided by the number of events.

Calculation of ΦFFCO 2 * surface fluxes
To calculate the surface fluxes, linear regression fits were made for each selected event. An example of a regression fit between FFCO 2 * and 222 Rn for a single event, which was observed from 13 August 2006 18:30 UTC to 7:00 UTC the next day, is shown in Fig. 4. 222 Rn increased by about 2 Bq m −3 during this period and FFCO 2 * by about 9 ppm. The correlation (assumed by the 222 Rn flux method) between FFCO 2 * and 222 Rn was significant: R = 0.9. For the regression slopes a correlation coefficient of R ≥ 0.7 is used as a threshold for further analysis. Hereby we aim to maintain regression slopes with a sufficient correlation and still have enough events to determine annually averages. In total 184 events were selected with 97 from sector 1 and 87 from sector 2. The surface flux FFCO 2 * for each event is calculated as follows (Levin, 1984;Thom et al., 1993;Schmidt et al., 1996) FFCO 2 * = Rn · FFCO 2 * Rn , where FFCO 2 * represents the surface flux of FFCO 2 * , Rn is the 222 Rn soil flux rate, FFCO 2 * represents the atmospheric mixing ratios of FFCO 2 * and Rn is the atmospheric concentration of 222 Radon. The overbars are the means for the spatial and temporal scales of influence. The 's represent enhancements from their background values. Hence, the second term in eq. 5 is the linear regression fit for each event. Basically, an atmospheric transport coefficient is determined from the ratio of the 222 Rn soil flux rate to the observed 222 Rn concentrations at a certain height. The mixing ratios of FFCO 2 * , if measured at the same height, can then be scaled to its surface flux ( FFCO 2 * ) assuming the atmospheric transport and dilution (e.g. by mixing with the free troposphere under unstable atmospheric conditions) is the same. As explained above, the radioactive noble gas 222 Rn is wherein τ max represents the maximum transit time for the specific observation period, for which we take the total length of each event, and λ is the decay constant (0.182 d −1 ) of 222 Rn. In this study, the applied correction for radioactive decay was on average <5%. For both sectors we used the 222 Rn soil flux ( Rn ) based on a 0.5 • × 0.5 • 222 Rn soil flux map for Europe (Szegvary et al., 2009) with a temporal resolution of 1 week. Szegvary et al. (2009) constructed this map, based on the correlation between 222 Rn and gamma-dose rate measurements based on observations in the year 2006. The annual mean 222 Rn soil flux was 0.29 atoms cm −2 s −1 for sector 1 and 0.35 atoms cm −2 s −1 in the case of sector 2. A very small seasonal cycle is present of about 0.01 atoms cm −2 s −1 with higher emissions in the summer and lower in the winter. The data was fitted with a weighted LOESS fit which we used for calculating the FFCO 2 * emissions (Fig. 5).

FFCO 2 surface emissions from modelled mixing ratios
To investigate our method's sensitivity towards emissions from sector 1 and 2 we performed the same exercise on modelled mixing ratios for the year 2007 from the regional atmospheric transport model REMO (Langmann, 2000). Mixing ratios of FFCO 2 and 222 Rn concentrations were simulated by REMO for station Lutjewad and from them FFCO 2 surface fluxes were calculated using the approach described in this paper for sector 1 and 2. The results are then compared to the a priori, annual mean fluxes for both sectors (Fig. 6). REMO has Fig. 5. 222 Rn soil flux rate (expressed in atom cm −2 s −1 ) for both sectors based on Szegvary et al. (2009). The data was fitted with a weighted LOESS to calculate the surface fluxes FFCO 2 * . a grid resolution of 0.5 • × 0.5 • in a rotated spherical coordinate system and was run on a semi-hemispheric domain covering the area north of 30N. For this study, REMO was fed with hourly FFCO 2 emissions based on 1 • × 1 • data available from the Emission Database for Global Atmospheric Research (EDGAR3.2FT2000) (Olivier et al., 2005). The emission data were extrapolated from the year 2000 to 2007 by C. Gerbig (personal communication, 2009) using BP statistics (downloadable from www.BP.com/statisticalreview) and seasonal to diurnal variations were included based on time profiles available in the EDGAR database. The 222 Rn soil flux was taken from the European map given by Szegvary et al. (2009) and has a temporal resolution of 1 week and a spatial resolution of 0.5 • × 0.5 • . This is the same source as used for the analysis of the observations (Fig. 5). Figure 7a shows the 222 Rn concentrations for the Lutjewad observations and the simulations using REMO/Szegvary for January-April 2007. The synoptic variations are very well captured by the model. Figure 7b shows the correlation between the 222 Rn concentrations from REMO/Szegvary and our observations for all selected events for the year 2007. The correlation between the two is good (R = 0.85) which suggests that the a priori 222 Rn soil flux map from Szegvary et al. (2009) is representative for the footprint of station Lutjewad. There is, however, a large scatter in the data, which may be related to the modelled atmospheric transport and/or the limited spatial resolution of the soil flux input data.   There are too few observations to determine the annual mean emissions per year. As expected based on the population density distribution, fluxes from sector 1 are generally lower then those from sector 2. The frequency distributions of the fluxes for all observations are shown in Figs. 9a, b (for the trajectories starting within section 1) and c (for those starting within sector 2). The results for the REMO based data are shown in Figs. 9d (all observations), 9e (sector 1) and f (sector 2). As explained in Section 3.2 we consider the emissions based on trajectories starting in sector 2 as representative for the emissions of the Netherlands whereas those starting in sector 1 represent the local emissions. The annual mean emissions representative for sector 2 itself are calculated by subtracting the local emissions for sector 1 from the total emissions, by taking 19 grid cells for the Netherlands, 11 grid cells for sector 2 and 8 grid cells for sector 1 (Fig. 6). For the total period of May 2006-June 2009 we find annual mean FFCO 2 * surface fluxes for the Netherlands (i.e.

Results
represented by emissions starting from sector 2) of: (4.7 ± 1.1) kt km −2 a −1 . For sector 1 we find: (3.3 ± 1.1) kt km −2 a −1 and this leads to (5.7 ± 1.1) kt km −2 a −1 for sector 2. For the REMO exercise we compare our results with the emission input into REMO for both sectors as indicated in Fig. 6. Using our approach we find annual mean FFCO 2 surface fluxes for the Netherlands of: (4.5 ± 0.4) kt km −2 a −1 , for sector 1 of (3.3 ± 0.6) kt km −2 a −1 and then (5.4 ± 0.6) kt km −2 a −1 for sector 2. These are in very good agreement with the input fluxes of: 2.9 kt km −2 a −1 for sector 1, 5.8 kt km −2 a −1 for sector 2 and 4.6 kt km −2 a −1 for the Netherlands in total. Hence: the 222 Rn flux method on the modelled concentrations of 222 Rn and FFCO 2 returns the input fluxes for both sectors. This is strong support for our chosen methodology: combining FFCO 2 * and 222 Rn observations with the presented analysis based on the selection of emission 'events' and back trajectory selection. Our methodology shows that using data from our station Lutjewad, emissions from both sectors and thus the Netherlands as a whole can be determined, provided a satisfactory number of 'events' with suitable back trajectories is available. All results, including our analysis with the modelled data from REMO, are summarized in Table 1.

Discussion
In this paper, we estimated the annual mean emissions of CO 2 from fossil fuels ( FFCO 2 * ) using the 222 Rn flux method on specific events which are selected based on their trajectories. In cases such as this one, it is common practice, but incorrect, to treat the data as shown in Fig. 9 as a normal or lognormal distribution with a single "true" value while the variability of the data is assumed to be caused by measurement and other random uncertainties. On the contrary, each measurement is the flux result for the specific area sampled by the air masses at that  4.5 ± 0.3 (4.5 ± 0.2) 4.7 ± 1.6 (kt km −2 a −1 ) specific event. Therefore, flux results for the individual events can vary widely and yet all of them should have their share in the result for the flux of the total area. The mean emission for a given area is thus represented best by the sum of all observations divided by the number of observations (i.e. the arithmetic mean), and its uncertainty is given by quadratic addition of the random errors of the individual observations. Our results regarding the national CH 4 and N 2 O emissions ( Van der Laan et al., 2009b), for which we used the same event-based 222 Rn technique, but for which we presented the results based on the median and a lognormal fit, will therefore be revised according to our approach presented here (i.e. using two sectors and the arithmetic mean) in a forthcoming paper together with an extended data set. The new results for the fluxes for the Netherlands so far are given in Van der Laan (2010). They are higher by about 50% compared to the lognormal approach and still confirm the national inventories.
Of course, the reliability of our method is crucially dependent on the completeness of coverage of the area, and thus on the total number and diversity of the events. Figure 3 is encouraging in the sense that it shows that our observations cover the footprint very well. Furthermore, the results we achieved from testing our method on the simulated data from the REMO model (i.e. calculating fluxes from individual events and analysing their trajectories), fed with EDGAR FFCO 2 emission data and the 222 Rn field of Szegvary et al. (2009) are crucial, and are encouraging: our method returns what was put into the model for each sector. This clearly supports our methodology.
Our results are potentially biased towards the warmer months since we have less observations in the winter months of 2007 and 2008 (Fig. 8). We expect this potential bias however not to be significant since the results for the last, in terms of observations more nicely distributed, year (August 2008-June 2009) would only be 11% lower.
The uncertainty for our observation-based emission estimate for the Netherlands as given in Section 4 was about 23% based on the random error, calculated as explained above. For the combined uncertainty we should also consider systematic errors, related to the conversion from 14 C to FFCO 2 * and the assumed 222 Rn soil emission rate. In this study the used 222 Rn soil flux was based on Szegvary et al. (2009). The uncertainty was estimated by the authors to be ±30% for the annually mean 222 Rn soil flux of the Netherlands. However, Figs. 7a and b show that the modelled 222 Rn concentrations agree very well with our observations, suggesting, at least for the annual means, that the 222 Rn soil flux is representative for the Netherlands. By using two sectors, we take into account part of the spatial variation and also a small seasonal cycle in the LOESS fit (Fig. 5) was included. On small time scales the uncertainty can potentially be large, which is mostly attributed to the fact that although the production of 222 Rn in the soils is constant its emission is influenced by atmospheric pressure, soil temperature and soil humidity. For example, the 222 Rn soil flux can vary by an order of magnitude due to, for example, rain or snow and proportionally affect the calculated fluxes. This probably explains, together with the uncertainty related to the modelled transport in REMO, the large scatter in Fig. 7b. However, we expect this effect to be a minor contribution in our results as we calculated the surface fluxes per single event of which the average duration is in the order of about 10 h. For our annual mean FFCO 2 * emissions we therefore estimate the uncertainty related to the 222 Rn soil flux at about 10%. The systematic uncertainty in our FFCO 2 observations is directly introduced into the dCO/FFCO 2 ratio for which we applied a weighted harmonic fit with a linear component to calculate FFCO 2 * . The fit is mostly determined by the last 2 yr of observations since there are far less observations in the first year, and the observations we do have are subjected to relatively large uncertainties. The uncertainty in the fit is estimated to be about 20% based on the uncertainties in the amplitude and the offset. However, this fit is based on 2-weekly integrated observations of CO and FFCO 2 whereas our events are typically 10 h. The ratio of dCO/FFCO 2 could be influenced by a diurnal variation related to the fact that different fossil fuel combustion processes (e.g. domestic heating and traffic) have very different CO/FFCO 2 ratios, and their relative contribution varies over the day (Levin and Karstens, 2007). Another potential source of uncertainty could be biomass burning or photochemical effects (Campbell et al., 2007). More research is needed to better estimate the temporal variation in the dCO/FFCO 2 ratio, for example by analysing events on a high-resolution basis for CO and FFCO 2 by using flask samples. We do not expect this uncertainty to contribute significantly to our annually averages, but our event selection procedures could potentially introduce a bias. Another source of uncertainty is the fact that the 222 Rn flux method is based on (vertical) atmospheric gradients which are observed mostly in the evenings and nights when the atmosphere is in general more stable (Fig. 10). Our method is therefore less suitable for estimating surface emissions in the afternoon when vertical mixing is more pronounced. Most of the day is, however, well covered and also the traffic peaks in the mornings and evenings are generally included in our data set. Figure 10 shows furthermore that there is no significant correlation between the height of the flux and the time of the events, which is probably because each event represents a single integrated value that usually includes emissions over several hours during day and night.
The sensitivity to the value of the correlation coefficient R as a threshold for the events was also tested: our results would be about 8% higher (but less representative due to the lower number of events) when choosing R ≥ 0.8 as a cut-off value, and 2% lower with R ≥ 0.6. Our results are therefore not significantly biased by our choice of R.
In total, we estimate the systematic error to be ±25%. Together with our random error of 23% this leads to a combined uncertainty for our annually mean results (Table 1) of ±35%. Fig. 10. Distribution of the selected events during the day. Fewer events are observed during 10:00-16:00 when vertical mixing is more pronounced.

Conclusions and outlook
The method presented here provides an important independent approach for validating inventories. We conclude that our method of 222 Rn flux approach leads to a reliable estimate (±35%) of the FFCO 2 fluxes for the Netherlands. For detecting year-to-year variations the uncertainty is lower (±23%). Our methodology is firmly backed-up by the simulated data treatment using REMO, for which our method reproduced the input emissions of EDGAR, and analysis of the distribution density of our observations towards the footprint. This implies that, using our method, a single monitoring station is in principle capable of determining the FFCO 2 flux for an area at least the size of the Netherlands (36 000 km 2 ). However, it is crucial that the target area is sufficiently covered by the trajectories of the collected events. Although in this sense the location of station Lutjewad is close to ideal, this requirement is ideally fulfilled with a network of monitoring stations to ensure that each grid on the total footprint is covered semi-continuously.
Based on 3 yr of observations from station Lutjewad for the period of May 2006-June 2009 we estimate net emissions for fossil fuel derived CO 2 of: (4.7 ± 1.6) kt km −2 a −1 for the Netherlands. Our result agrees very well with the Dutch inventory of (4.5 ± 0.2) kt km −2 a −1 (average value from 2006 to 2008) but more research is needed to reduce the uncertainties.
As a side result we found that the 222 Rn soil flux estimates of Szegvary et al. (2009) are representative for the Netherlands. Still, for future studies we suggest more research to further constrain the 222 Rn soil flux. If future research would result in a more reliable soil 222 Rn emanation rate, our values can be adjusted in a simple way as they are directly proportional to this emission rate. We encourage that our method is applied using data from other stations. We suggest that the sensitivity towards specific regions (i.e. how each grid is covered) is tested with a higher resolution forward transport model in combination with a high-resolution emission inventory (e.g. the newly available 0.1 • × 0.1 • EDGAR4.0 database). For each individual event the surface fluxes could then be calculated based on its individual trajectory.