Estimates of European emissions of methyl chloroform using a Bayesian inversion method

. Methyl chloroform (MCF) is a man-made chlorinated solvent contributing to the destruction of stratospheric ozone and is controlled under the “Montreal Protocol on Substances that Deplete the Ozone Layer” and its amend-ments, which called for its phase-out in 1996 in developed countries and 2015 in developing countries. Long-term, high-frequency observations of MCF at three European sites show a constant decline in the background mixing ratios of MCF. However, we observe mixing ratio enhancements of MCF in pollution we long-term

For this reason MCF was included in the 1990 amendment to the Montreal Protocol, requiring it to be completely phased out in 2005 and 2015 in non-Article 5 parties and Article 5 parties, respectively. The phase-out in non-Article 5 parties was significantly accelerated (from 2005 to 1996) by the 1992 amendment to the Montreal Protocol (UNEP, 2009). Due to the relatively short lifetime of MCF and its past use primarily in fast-release applications, the implementation of the Montreal Protocol allowed, starting from 1998, a nearconstant exponential-rate decline in the global background concentrations of MCF (Rigby et al., 2013;Montzka et al., 2011a). Currently, the global average background mixing ratio of MCF is around 5 ppt, which can be compared to the maximum value of ca. 130 ppt reached in the mid-1990s (Montzka et al., 2011b). However, latitudinal and vertical mixing ratio gradients, as reported by Montzka et al. (2000) on the basis of analysis of flask samples collected across the globe, still persist, consistent with continued northern hemispheric emissions or an asymmetry in MCF loss in the Northern vs. Southern Hemisphere (data update to Montzka et al., 2000 available at ftp://ftp.cmdl.noaa.gov/hats/solvents/ CH3CCl3/flasks/GCMS/CH3CCL3_GCMS_flask.txt).
Based on our understanding that OH is the main removal mechanism for MCF in the troposphere, long-term global measurements of MCF, combined with emission estimates, have been used to infer both hemispheric and global average mole fractions and trends, as well as variability of the OH radical on annual timescales (see, for example, Montzka et al., 2000Montzka et al., , 2011aPrinn et al., 2001Prinn et al., , 2005. This is otherwise difficult because no globally representative measurements of OH exist, although OH plays a key role in atmospheric chemistry, due to its capability of removing many other radiatively important trace gases and atmospheric pollutants (Montzka et al., 2011a).
A correct quantification of MCF emissions is crucial for evaluating hemispheric and global mean OH concentrations. In recent years, this has stimulated a scientific debate concerning the actual extent of emissions, in particular on a regional (European) scale. According to Krol et al. (2003), who used data from a short-term tropospheric measurement campaign, annual European emissions in 2000 were greater than 20 Gg, whereas Reimann et al. (2005), based on 4 years (2000)(2001)(2002)(2003)(2004) of high-frequency data from two European stations (Mace Head, Ireland, and Jungfraujoch, Switzerland), have derived considerably lower annual emissions, ranging from 0.3 to 3.4 Gg yr −1 .
Long-term, high-frequency observations of MCF are regularly conducted at four sites in Europe that are part of the former SOGE (System for Observations of halogenated Greenhouse gases in Europe) network: Mace Head (MHD, Ireland), Jungfraujoch (JFJ, Switzerland), Zeppelin mountain near Ny-Ålesund (ZEP, Norway) and Mt. Cimone (CMN, Italy). The network was established in 2001, though measurements at MHD had started many years earlier , with its aim being the creation of a European in-frastructure for the evaluation of long-term trends and European emissions of halogenated gases. In fact, long-term, high-frequency observations combined with inverse modelling (top-down approaches) have proved to be a powerful and important tool for the quantification of emissions and the verification of bottom-up inventories for many trace gases (Nisbet and Weiss, 2010).
The time series of MCF mixing ratios from the European stations show a constant decline that is consistent with the global trends in atmospheric mixing ratio. However, the persisting occurrence of non-negligible mixing ratio enhancements during pollution episodes, especially at CMN and JFJ, suggests unexpectedly high ongoing emissions somewhere in Europe. In order to identify the source region and estimate the magnitude of such emissions, we have used a Bayesian inversion method developed by Seibert (2000Seibert ( , 2001 and improved by Eckhardt et al. (2008) and Stohl et al. ( , 2010 as well as a point source analysis method developed by Keller et al. (2011), in combination with long-term, highfrequency observations at the three European stations CMN, JFJ and MHD. The estimates provided by this analysis are relevant for constraining the atmospheric budget of MCF on a regional scale, improving the accuracy of MCF emissions quantification on a global scale and, ultimately, improving our understanding of atmospheric hydroxyl.

Measurement sites
The measurement sites CMN (2165 m above sea level, a.s.l.) and JFJ (3580 m a.s.l.) are continental stations located in the northern Apennines in Italy and in the central Swiss Alps, respectively. MHD is a marine station located on the western coast of Ireland (Grant et al., 2010). The three stations are remote measurement sites episodically influenced by emissions from the European continent (Fischer et al., 2003;Henne et al., 2010) and are connected to the international measurement network AGAGE (Advanced Global Atmospheric Gases Experiment). Furthermore, they are classified as WMO GAW (World Meteorological Organisation Global Atmosphere Watch) global stations.
Time series from ZEP have not been included in the analysis because tests performed on 1 year of data (2008) showed that the inclusion of this remote station, where no MCF enhancements are observed, did not substantially alter the emission magnitudes we derive here (derived fluxes differed by < 2 % with the inclusion of ZEP).

Analytical method
At CMN, bihourly observations of MCF in ambient air are conducted via gas chromatography-mass spectrometry (GC-MS) preceded by online sample enrichment using adsorbent material. An Agilent 6850-5975 GC-MS system has Atmos. Chem. Phys., 14, 9755-9770, 2014 www.atmos-chem-phys.net/14/9755/2014/ been equipped with an auto-sampling/pre-concentration device (Markes International, UNITY2-Air Server2) to enrich the halocarbons on a focussing trap filled with four different adsorbing materials. The details of the analytical procedure are reported in Maione et al. (2013). The MS detector is operating in selective ion mode (SIM), and m/z values of 97 and 61 are selected for detection and quantification of MCF. Typical repeatability evaluated from the repeated working standard measurements for MCF is 1.8 and 1.3 % (relative standard deviation, RSD, over 1 year), referring to the period before and after 2008, respectively, i.e. before and after the analytical instrumentation (mass spectrometer and pre-concentration unit) was upgraded so as to improve the quality of data, with a clear decrease in instrumental noise. The relatively lower precision in the first half of the data set does not affect data reliability. Limit of detection, LOD (signal-to-noise ratio, S / N > 3), and limit of quantitation, LOQ (S / N > 10), have been indirectly estimated on the chromatogram of the working standard run with typical background concentrations. The obtained values are the following: LOD 0.40 and 0.30 ppt and LOQ 1.4 and 0.9 ppt, before and after 2008, respectively. The JFJ and MHD MCF time series used in this paper are based on an ADS (adsorption-desorption system; Simmonds et al., 1995) and MEDUSA (Miller et al., 2008) measurements. At JFJ, measurements started in January 2000 using the ADS, and in 2008 this sampling device was replaced with a MEDUSA system; at MHD, measurements started in 1998 using the ADS, and in 2004 it was replaced with a MEDUSA system. Typical repeatability for the AGAGE MEDUSA systems is well below 1.5 % for both stations (RSD over 1 year).

Calibration
Ambient air samples are analysed every 2 h and are bracketed by working standard analyses in order to detect and correct for short-term instrumental drift. These working standards are air samples which are pumped during relatively cleanbackground conditions into 35 L electro-polished stainless steel canisters (Essex Cryogenics, Missouri, USA) using a modified oil-free compressor (SA-3 or SA-6 Rix, California, USA) up to ∼ 60 bar. The tanks are humidified with purified water during the pumping process in order to improve the stability of the compounds (Miller et al, 2008) and to ensure a close similarity in composition between the ambient air samples and the reference standard. This with the aim of minimising analytical artefacts and the nonlinearity of the method. The working standards are calibrated at least monthly against a tertiary standard prepared with the same procedure at the baseline station of MHD. For JFJ and MHD, transfer standards from SIO (Scripps Institution for Oceanography, La Jolla, CA, USA) are used for this purpose. The calibration scale for MCF used here is the SIO-2005 scale. The GC-MS systems at the three stations are fully automated via the Linux-based chromatography soft-ware (GCWerks, gcwerks.com) adopted by the AGAGE programme.

Determination of concentrations in the background atmosphere
A careful evaluation of MCF mixing ratios in the background atmosphere reaching our measurement locations is crucial not only for estimating atmospheric trends and, consequently, annual growth rates but also for emission evaluation, because back attribution techniques used for assessing emissions are based on the clear identification and quantification of mixing ratio enhancements above background values.
In the case of a station like CMN, surrounded by complex topography and emission fields, the determination of the baseline is not trivial. Here we apply a statistical method based on a two-step procedure (Giostra et al., 2011) that we have developed specifically for CMN and also applied to mixing ratio histories acquired at JFJ and MHD. The first step consists of detrending the measurement data record using an appropriate time interval. The second step is aimed at estimating the uncertainty in the determination of the baseline, which includes instrumental error and natural background variability, assuming that such errors follow a Gaussian distribution. The overall observed probability distribution function (PDF) can be decomposed into the sum of a Gaussian and a gamma distribution with the Gaussian distribution corresponding to the well-mixed background atmosphere state and the gamma distribution corresponding to a non-well-mixed state, i.e. data containing recent emission inputs that will be used to derive information about regional emission rates and their spatial distribution. If the number of available data points is large enough, then the decomposition of the global PDF for the observations as the sum of a Gaussian plus a gamma becomes stable and reliable. Data belonging to the obtained Gaussian distribution are regarded as baseline data.

Dispersion modelling and Bayesian inversion
The inversion procedure is based on 20-day backward simulations with the Lagrangian particle dispersion model FLEX-PART (Stohl et al., 1998(Stohl et al., , 2005Seibert and Frank, 2004).
FLEXPART is a stochastic model with detailed treatment of turbulence and convection and uses meteorological analyses from the European Centre for Medium-Range Weather Forecasts (ECMWF). In this study we used the ECMWF analyses at 1 • × 1 • resolution for the period 2002-2012, over the domain reported in Fig. 1. In addition, over a 2-year period (2008-2009), we used nested meteorological data with a resolution of 0.25 • × 0.25 • in the European domain (from 12 • W to 28 • E and from 35 to 65 • N), called ECMWF_nest. FLEXPART was run backward in time from the measurement stations at 3-hourly intervals, using 40 000 particles for each backward run. The FLEXPART output is an emission sensitivity, also called source receptor relationship (SRR). The SRR in a particular grid cell, expressed in units of s kg −1 , is proportional to the particle residence time in that cell and measures the simulated mixing ratio at the receptor that a source of unit strength (1 kg s −1 ) in the cell would produce for a given air sample. Multiplying the footprint emission sensitivity (i.e. the emission sensitivity in the lowest model layer) by the emission flux taken from an appropriate emission inventory gives the simulated mixing ratio at the receptor, which can be compared with a coincident measurement. Figure 1 shows the average footprint emission sensitivity for the three stations (CMN, JFJ, MHD) for the period January 2008-December 2009. Hereinafter, we will refer to the area with sensitivity > 2 ps kg −1 as the study domain.
The FLEXPART output can be ingested directly by the inversion algorithm based on the inversion method developed by Seibert (2000Seibert ( , 2001 and improved by Eckhardt et al. (2008) to allow for (i) an a priori emission estimate for the unknown magnitude and location of sources, (ii) a Bayesian formulation considering uncertainties for the a priori emissions and the observations, and (iii) an iterative algorithm for ensuring a solution with only positive values.
A further improvement was introduced by Stohl et al. ( , 2010, which considers a baseline in the observations that is adjusted as part of the inversion process, as well as more detailed quantification of errors. The model setup and the inversion method used for this study are the same as in Stohl et al. ( , 2010, where additional mathematical details can be found.
The basic idea is to find an a posteriori emission distribution leading to the best fit between the measurements and the model results while keeping the solution within the given error bounds of the a priori emissions. The "best" agreement is measured as the sum of the squared errors, inversely weighted with the uncertainty variances. The method used also identifies "outliers" in the model-simulated mixing ratios and assigns them large uncertainties to prevent the solution being strongly influenced by large measurement and/or model errors .
The a priori MCF emissions used for this study are based on the E-PRTR (European Pollutant Release and Transfer Register) inventory -a Europe-wide register that provides environmental data from industrial facilities in EU member, states, as well as in Iceland, Liechtenstein, Norway, Serbia and Switzerland and reports MCF atmospheric emissions higher than 100 kg yr −1 -from 2007 to 2011. For those years in which emission data are not available (i.e. 2002-2006 and 2012), we used the average emission values over 2007-2011. An alternative homogeneous emission field was used with aim to test model performance. Results are shown in Appendix B. The size of the inversion problem is defined by the number of grid cells for which emissions shall be determined. In order to reduce the number of unknowns we used a variable-resolution emission grid, with grid sizes ranging from 1 • × 1 • lat-long to 2 • × 2 • lat-long (Spain, Portugal and Norway). The resolution was controlled by the product of a priori emission flux and average emission sensitivity, as described by , who used 1 • × 1 • as the finest resolution. SRR values are high in the vicinity of the observation sites and they decrease with the distance from the sites (see Fig. 1), since emissions at large distances from the measurement locations cannot be resolved at high spatial resolution. Since MCF emissions predominantly occur over continents, in the inversion we excluded boxes that are covered by water or ice by more than 99 %, a higher value with respect to the 95 % value used in a recent study (Keller et al., 2012). The number of grid cells used for the annual inversions is 4400.

Point source analysis
The point source analysis (PSA) is based on the approach developed by Keller et al. (2011) called "(time variable) point source analysis" that attributes excess concentrations, not explained by the a priori emission field, to a source area whose location and extension is assumed to be known and that appears to be responsible of the enhancements. Given the emission sensitivity in the source area as obtained from the dispersion model, the source region emission magnitudes can be directly determined for each individual excess concentration event. This has the advantage that time-varying emissions can be retrieved, whereas the emissions are assumed to be constant during the course of a year for the inverse modelling. In the case of the point source analysis, Keller (2011) defines the sensitivity as the residence time of the particles below 100 m within a 1 • × 1 • cell, above the model ground, divided by air density, hereinafter defined as SRR v (SSR v : SRR volume; volume = 1 • × 1 • × 100 m). In this study, the volume is 0.5 • × 0.5 • × 100 m. When the emission sensitivity in the source region is low, spuriously high emission values can occur with this method. Several tests were therefore performed in order to assess, for the identified source area, a SRR v threshold above which these noise problems are largely avoided.
The PSA method takes each individual measurement above the baseline, not explained by a priori emission field, and determines the emission that is needed in a predefined source region to reproduce this measurement exactly. Due to inaccuracies in transport and other numerical errors, this normally results in a noisy emission time series that, upon averaging in time, provides an estimate of the emission in the predefined region. The method is particularly suited for cases when source region emissions are variable and also when it is not certain that emissions are constant.

Time series
The statistical filter as in Giostra et al. (2011) for the identification of baseline mole fractions has been applied to the time series data available at the three stations. The resulting time series are reported in Fig. 2, where the black dots represent baseline data and the red dots represent air samples with mixing ratios above the baseline. The monthly mean baseline mole fractions of MCF at the three stations have been used for the evaluation of the trends and the changes in the trend, using the empirical model as reported by Simmonds et al. (2004), whose results are reported in Table 1. MCF average baseline mole fractions show a significant and consistent decrease over the study period at all the three stations as a consequence of a strong decrease in global emissions following the implementation of the Montreal Protocol in a manner consistent with data obtained at other sites around the globe (Montzka et al., 2011a). However, the percentage of enhancements above the baseline, as well as their intensity, is noticeably lower at MHD compared to the other two stations, suggesting the persisting occurrence of recent emissions, which continue to influence the continental stations more strongly than MHD. It is thus likely that the MCF source region is located closer to and more upstream of CMN and JFJ than to MHD.

Bayesian inversion setup
In order to identify MCF source regions and quantify the magnitude of emissions responsible for the observed mixing ratio enhancements above background levels, the Bayesian inversion method was applied to the MCF data from the three sites. The a priori emission field was based on the atmospheric emissions from the E-PRTR data set, reported in Fig. 3. The E-PRTR database also includes MCF emissions to the soil and water, which have not been included in the a priori field.
Initially the simulations have been performed using observations from January to December 2008, for which nested meteorological data with a resolution of 0.25 • × 0.25 • are available.
The uncertainties of the emissions, σ x , need to be specified for every grid cell. As there is no information about uncertainties, we used the uncertainty in the matrix diagonal elements as defined in , i.e. for inversion box j , σ j x = max p · x j , 2 · p · x surf , with p being a properly chosen scaling factor; x j the a priori emission flux in the inversion box j ; and x surf the global emission value, as estimated by Rigby et al. (2013), homogeneously distributed in the grid cells corresponding to land areas. This σ formulation allows us to assign a high uncertainty also to those boxes with zero a priori emissions. We tested p values ranging from 50 to 500 % of the prior emission estimate, balancing between (i) enough flexibility in emissions to allow for adjustments that better fit the observations and (ii) not too high flexibility that might lead to over-fitting of the observations and to noisy and unrealistic emissions. The a posteriori flux assigned to the inverted box is distributed according to the population density.
However, the improvement of the measurement-model correlation was largest in the p value range from 0.5 to 1, with much smaller improvements for further increases of p. Furthermore, the retrieved emission field subjectively appeared too noisy for the highest p values. We therefore set p = 1 as a compromise, as this produced a reasonably low noise level in the a posteriori emissions while still obtaining a good correlation between observed and a posteriori modelled data. Sensitivity tests performed in order to investigate how the a priori emission intensity and the station network geometry affect the a posteriori emission field are reported in Appendixes A1 and A2.

Inversion results
The E-PRTR emission map and the resulting a posteriori emission field for 2008 are shown in Fig. 4, proving that the Bayesian inversion is able to (i) to confirm the localisation of the majority of the sources declared in the E-PRTR inventory as atmospheric emissions; (ii) localise sources not included in the a priori field, but included in the E-PRTR inventory as release to soil/water, e.g. in Norway; (iii) identify additional emissions sources not reported by E-PRTR, e.g. in northern Italy, where a large waste water treatment plant is present; and (iv) confirm that the strongest sources of MCF in Europe are located in southeastern France, hereafter named SEF. We define the SEF as the area including the two French sources declared in the E-PRTR inventory and that shows an emission more than 1000 % higher than adjacent cells. It is worth nothing that SEF includes the two strongest MCF sources in Eu-rope, as declared by E-PRTR (landfill and waste disposal or recovery of hazard waste and halogen chemical plant). Due to the importance of the SEF area within the study domain, we performed the inversions for the whole study period, isolating SEF as a separate region.
As shown by the sensitivity tests (see Appendix A), the 1 • × 1 • resolution and the nested 0.25 • × 0.25 • meteorological data produce comparable emission estimates. Therefore, we extend the proposed analysis to the whole period, January 2002 to December 2012, for which only the coarser resolution ECMWF data were available. The detailed 2002-2012 a posteriori emission estimates, and the related uncertainties (calculated as described in Appendix A1, Sensitivity tests) for different European areas are given in Table 2, where the a priori emissions are reported as well. These results confirm the importance of the SEF area, but they suggest an MCF emission roughly 5 times higher than the E-PRTR inventory estimate and up to 1 / 10 (in 2009) of the total semi-hemispheric (30-90 • N) emissions for the same years reported by Rigby et al. (2013). All other sources in Europe together contribute only about 50 % of the total European emissions.
In Fig. 5 (top panel), the trends in emission estimates from SEF, from different groups of countries, and from the whole European study domain are reported. Although the study domain and the SEF area both show a decrease in emissions over the whole study period, the decrease is less rapid for SEF and the estimated emissions from the SEF region are remarkably large. The extent of such estimates can be further appreciated by comparing our estimates from SEF (red curve) with the E-PRTR data for SEF used for the a priori (light blue) and the global (green) and semi-hemispheric (30-90 • N) (purple) emissions provided by Rigby et al. (2013), reported in the bottom panel in Fig. 5. The fraction of global emissions coming from the SEF region ranges from 2.6 % in 2003 to 10.3 % in 2009, with an average of 6 % over the whole period (black line in the bottom panel of Fig. 5). There also appears to be an upward trend in the fraction of global MCF emissions coming from SEF, indicating that MCF emissions in SEF decrease more slowly than in the rest of the world. These results confirm the relevance of the SEF region that accounts for the majority of the emissions in the study domain. The European study domain emissions are a significant fraction of the total semi-hemisphere emissions, also including Article 5 countries (Rigby et al., 2013), ranging from 9.8 % in 2004 to 33.7 % in 2011, of which on average 50 % are from the SEF region. The plot also shows that, as discussed in Appendix A, the a posteriori emission estimate is not severely impacted by the substantially different emission magnitude reported during 2010 for SEF in the E-PRTR inventory.

Meteorological filter
Since SEF has been identified as the strongest MCF source region in the study domain, a further analysis has been performed in order to highlight how emissions from this region are recorded at the receptors or sampling locations. Using the FLEXPART model output, the measurement data were categorised into cases with zero emission sensitivity in the SEF region (SRR SEF ) and cases with SRR SEF > 0.
The data series reported in Fig. 6 show that, at CMN and JFJ, nearly all mixing ratio enhancements are associated with SRR SEF > 0, and most data on the high side of the background distribution also have SRR SEF > 0. Even at MHD, which is a considerable distance away from SEF, all samples with SRR SEF > 0 are on the high side of the baseline distribution. This suggests that the SEF region is the dominant source contributing to mixing ratio enhancements at CMN and JFJ. Results from MHD are not inconsistent with this conclusion. However, since not all data at the high end of the baseline distribution are explained by non-zero SRR SEF , these data are likely to be a contribution of weaker sources within the study domain.
The detrended data set is reported in Fig. 7 in the form of probability density function (PDF) distributions, where the data pertaining to SRR SEF = 0 follow a Gaussian distribution, corresponding to a well-mixed state (baseline), and the data pertaining to SRR SEF > 0 (having non-zero sensitivity to the SEF region) follow a gamma distribution with a long tail towards high values, corresponding to the influence of recent emissions. The sigma (σ ) value of the distribution decreases with the increase in the extent of mixing, converging towards the instrumental error in the limit of "perfect" mixing. Here, σ values are 0.16 ppt at CMN, 0.17 ppt at JFJ and 0.15 ppt at MHD. The instrumental error at these different sites is 0.13, 0.15 and 0.15 ppt, respectively, during the period 2008-2009. The slightly higher σ values at CMN and at JFJ confirm that a spatial gradient for MCF is still occurring, confirming the occurrence of fresh emissions in continental Europe. For more details on the influence of the spatial gradient on the PDF, see Giostra et al. (2011).

Source quantification: point source analysis
In order to verify the emission magnitudes derived with the Bayesian inversion, we used an alternative method, called point source analysis, introduced by Keller et al. (2011), which is suitable once an emission point (or an emission area, in our case the SEF area) is identified. The advantage of this method, hereinafter called PSA, is that it does not rely on the assumption of constant emissions over a certain period. Another interesting feature of the application of the PSA is that independent estimates of emissions for the area can be derived with this analysis for each of the observing stations. The analysis has been conducted for the 2-year period (January 2008-December 2009) when the high-resolution meteorological data were available.
For the PSA, data have been selected when the sampled air was influenced by the SEF area. We set a threshold when the SRR v in that area was above 1500 s m 3 kg −1 . During those periods, excess concentrations that cannot be explained by the a priori emission field are assumed to originate exclusively from the SEF area. For every measurement point with a sensitivity above 1500 s m 3 kg −1 , a SEF emission flux is derived that will bring the modelled mixing ratio (given the simulated SRR v values) into perfect agreement with the measurements. The use of a minimum threshold is necessary to avoid the greater uncertainties associated with sampling events characterised by very small sensitivities in the SEF region.
The SRR v threshold (1500 s m 3 kg −1 ) is identified as the value for which the annual emission flux converges towards a Atmos. Chem. Phys., 14, 9755-9770, 2014 www.atmos-chem-phys.net/14/9755/2014/  limit value. The plots in Fig. 8 (top panels) report the average annual emission fluxes together the 95 % confidence interval of the annual average as obtained by 3-hourly estimates of the three stations. The bottom panels report the number of measurements per year above the identified threshold. The emission estimates from the SEF region are 0.5 ± 0.1 and 0.5 ± 0.2 Gg yr −1 in 2008 and 2009, respectively. The obtained estimates using the PSA method are almost a factor of 2 larger than those obtained from the inverse modelling, as reported in Table 2. This is expected, as the inversion results are bound towards a priori emissions that are clearly too low for the SEF area. This leads to a low bias also in the obtained a posteriori emissions.

Additional analysis
Our results suggest that the SEF area is responsible for a large fraction of semi-hemispheric MCF emissions. To further as- sess this hypothesis, we consider additional measurements during a sampling campaign carried out at four sites located in the vicinity of the Marseille urban area (well inside the SEF region) in June and July 2001. The obtained time series is reported in Fig. 9, along with the MHD baseline time series for the same period (grey dots). Comparing MCF mixing ratios measured in the Marseille area with the MHD record, it is possible to appreciate the magnitude of the enhancements above the Northern Hemisphere baseline (represented by the MHD record). MCF mixing ratios show a large variability with about 6 % of the data above 50 ppt and about 2 % of data above 100 ppt. Such large anomalies indicate the presence of a strong MCF source. Indeed, MCF data from several urban areas taken from Barletta et al. (2006) and reported in Table 3 show that the mean MCF mixing ratios in Marseille are among the highest during these years and, particularly, that the standard deviation of the data is by far the largest of all the sampled cities. This must be attributed to fresh emissions from nearby sources. The number of available data in this short localised sampling campaign is not large enough to allow for a precise localisation of the MCF source. However, we produced numerical simulations with a simple dispersion model, which was based on local wind data at Marseille. We used a source strength that was extrapolated to 2001 from our 2002-2012 emission estimates from the SEF area (i.e. 1.1 Gg yr −1 ) and tested various source locations in the Marseille area. The results obtained, albeit qualitative, confirm that the Table 3. MCF levels measured in various cities worldwide (adapted from Barletta et al., 2006).

Sampling site
Mixing ratio (ppt) 1σ standard deviation (ppt) Study period enhancements measured in the Marseille area in 2001 are compatible with a strong source of the magnitude obtained from our PSA located within some few tens of kilometres from the sampling sites.

Discussion and conclusions
MCF measurements carried out at three European sites (CMN, JFJ and MHD) from January 2002 through December 2012 show a number of regular and surprisingly persistent enhancements above the baseline levels. In order to identify and quantify potential European source regions for MCF, we set up a Bayesian inversion methodology initially tested over a 2-year period (2008-2009) for which high-resolution (0.25 • × 0.25 • ) meteorological data were available, which was then subsequently extended to the full 11-year record with lower resolution meteorological information.
The a priori emission field was based on the E-PRTR inventory, but the regions identified as strongest sources were also apparent with a less accurate a priori field (homogeneous; see Appendix B). Most of the sources declared in the E-PRTR inventory are landfills and disposal or recovery sites of hazardous waste. However, there is a strong inhomogeneity among the number of sources declared by the European countries reporting to E-PRTR, with the UK declaring more than 30 sources out of a total of 54, while several countries (e.g. Italy) do not report any emission. Since most of the sources declared in the UK are not linked to industrial processes but instead to population, it seems that the information from other countries is lacking. The methodology used allowed us to localise not only the source regions included in the a priori field but also additional sources such as sources in northern Italy, and it also highlights the model performance throughout the study domain. Furthermore, the method allowed us to identify sites where MCF is released to soil and water but the loss into the atmosphere is not given. These particular sites have therefore not been included in our a priori emission field. This could suggest that these releases to water and soil ultimately, at least partly, also end up as emissions into the atmosphere. However, since these source regions are far away from the receptors and are probably not particularly intense, the method is not able to give a quantitative estimation of the emissions.
Our results confirm that the strongest sources of MCF in Europe covered by our study domain are located in southeastern France (SEF), where landfill and waste disposal and halogen chemical plants are present. The quantification of low-intensity emissions, outside the SEF region, includes an uncertainty of around 100 %, and so these results cannot be used to suggest that the declared emissions are incorrect. However, the main European source region, the SEF area, is well defined with a 30 % uncertainty, and these estimated emissions are 5 times larger than the value declared in the E-PRTR inventory.
Such estimate was validated by a point source analysis applied to the same 2-year record, whose results provide emission estimates higher than the Bayesian inversion, confirming that the E-PRTR database underestimates emissions from the SEF region.
The estimated emissions, even though showing a strong decay over the whole study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)  Following Stohl et al. ( , 2010 various sensitivity tests have been performed throughout the entire study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) in order to investigate how the a priori emission intensity and the station network geometry affect the a posteriori emission field. In addition, for the years 2008-2009, for which alternative meteorological data are available, we investigated the influence of the meteorological data resolution. These tests provided a set of estimates, whose average is our best estimate, which is affected by an uncertainty corresponding to the maximum error. The maximum error is defined as the semi-difference between the maximum and minimum value of the a posteriori emission fluxes. In addition, we derived the percentage ratio (R p ) between the maximum error and average emission values.

A1 A priori emission field modulation
We scaled the a priori emission field with four different values: 0.5, 1, 1.5 and 2. Despite the differences in the emissions by a factor of 4, the a posteriori SEF emissions remained nearly constant, with an R p between 0.6 and 16 %. Emission estimates from areas other than SEF, however, show higher variations, with R p ranging from 20 % to more than 100 %. This is due to very low emission in those areas that are not well constrained using this approach. The histogram in Fig. A1 shows the ratio between the maximum difference among the a posteriori emissions obtained for the four error test cases and their average value for the year 2008. For the entire investigated period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), the average R p is only 9 % for the SEF area, suggesting a high reliability of the emission estimate, whereas the average values are much larger for areas other than SEF, reaching almost 100 %, thus indicating uncertainties of the same order of the retrieved emissions. The 2008 and 2009 R p values for the SEF region are reported in Table A1.

A2 Station network geometry
The effect of the station network geometry on the inversion results has been tested by removing one station at a time, rerunning the inversion and comparing it to the reference inversion. In this case the SEF region shows an averaged (2002-2012) R p of 18 %, suggesting that the inversion results, with our testing, are more affected by the geometry of the station network rather than the a priori field modulation. For the other regions the averaged relative differences is around 40 %. R p values for all regions and for 2008 are shown in Fig. A2.

A3 Meteorological data resolution
Finally, the relative differences between a posteriori emissions were calculated for the inversions based on FLEX-PART model runs using meteorological input data from ECMWF with different spatial resolutions: 1 • × 1 • and 0.25 • × 0.25 • . The R p values obtained with two different wind field spatial resolution for 2008 and 2009 relative to the SEF area were 3 and 2 %, respectively (Table A1). The overall uncertainty associated with the estimated emission values are calculated from the overall relative uncertainty derived from above described sensitivity tests. The errors associated Table A2. Single station parameters. Mean; SD, standard deviation; N , number of observation; E a , rms error priori; E b , rms error posteriori; 1−E a /E b , relative error reduction; E n b , a posteriori error normalised with the standard deviation of the observed concentration minus baseline; r 2 a , squared Pearson correlation coefficients between the observations and the total a priori; r 2 b , squared Pearson correlation coefficients between the observations and the total a posteriori; r 2 ba , squared Pearson correlation coefficients between the observations and the a priori baseline; r 2 bb , squared Pearson correlation coefficients between the observations and the a posteriori baseline; r 2 ea , squared Pearson correlation coefficients between the observation minus the a priori baseline and the modelled a priori; r 2 eb , squared Pearson correlation coefficients between the observation minus the a posteriori baseline and the modelled a posteriori.

A4 Single station parameters
Another way to evaluate the inversion performances is to analyse the time series statistical parameters for each station. As in , the station-specific error statistics were evaluated by comparing the a posteriori and a priori errors at different stations. The relative error reduction 1−E b /E a values (see Table A2) for CMN and JFJ was 0.15, and for MHD was 0.30, showing that the error reduction for the mountain stations is significantly lower than for MHD, being related to transport episodes associated with localscale wind systems that cannot be captured by the model. During such episodes, model errors cannot be reduced by improved emission data, yielding lower overall error reductions. As shown in , the two mountain stations are affected by incurable errors due to the complex topography, while the dispersion model shows the best performance in flat areas, where the meteorology is well described by the ECMWF data.
A substantial part of the MCF signal observed at the stations is explained by the variability and trend of the baseline, expressed as the squared Pearson correlation coefficient, r 2 ba , between the a priori baseline and the observed concentration (Table A2). As the a priori and a posteriori baselines are quite similar, the statistical results for both are nearly identical. As shown in , r 2 ba is higher for remote stations where events with transport from source regions on the timescale of 20 days are rare (e.g. r 2 ba = 0.91 for MHD) and intermediate at stations not too far from source regions (e.g. r 2 ba = 0.74 and r 2 ba = 0.55 for JFJ CMN, respectively) where the short-term variability is large.
The variability of values above the baseline reflects the occurrence of pollution events affecting the measurement sites. The model capability to capture these events is described by the correlation analysis of the polluted events with the simulated emission contributions from the last 20 days, using either the a priori (r 2 ea ) or the a posteriori model results (r 2 eb ) (see Table A2).
The three stations used for the inversion present rather low r 2 ea and r 2 eb values -the two mountain stations because of the incurable errors reported above, and MHD because it is far away from the main source region, located in SEF, and because measuring values that are typical also for remote stations such as Samoa (see, for example, ). However, for all stations, r 2 eb values are increased compared to r 2 ea values, suggesting that the a posteriori emission field is closer to the real emission field than the a priori emission field.
Atmos. Chem. Phys., 14, 9755-9770, 2014 www.atmos-chem-phys.net/14/9755/2014/ Appendix B Figure B1. A posteriori study domain emissions in pgm −2 s −1 of MCF based on measurements at three European sites and using the BI-1 inversion. Reference period: January-December 2008. White crosses over black dots indicate the locations of the measurements sites. Figure B2. As in Fig. B1 but using the BI-2 inversion.
In order to investigate the algorithm performance, we run the inversion without any information on source distribution in the a priori emission field. A detailed description of this test is reported in Graziosi (2013). An a priori homogeneous emission field was chosen based on the semi-hemispheric (30-90 • N) emission values reported in the study by Rigby et al. (2013), which provides global and semi-hemispheric emission data for the period 1951-2013 derived from atmospheric concentration data from the NOAA and AGAGE networks. The Bayesian inversion method was applied to the MCF time series available at the three sites. The simulations were performed using data from January 2008 to December 2008, when nested meteorological data with a resolution of 0.25 • × 0.25 • were available. The uncertainty of the emissions, σ x specified for every grid cell, σ j x = px j , is set using p = 1, corresponding to an uncertainty of 100 %. We refer to this setting as Bayesian inversion 1 (BI-1).
The resulting a posteriori emission field for 2008 is reported in Fig. B1. A clear hot spot is highlighted in SEF, in a region that is some 250 km × 250 km large, corresponding to the main source declared on E-PRTR inventory. Additional small and intense source regions are present in northern Norway, the southern UK and Benelux, as well as in northern Italy. All these sources have a small intensity compared with the SEF area. Most of them are reported in the E-PRTR inventory, while sources in northern Italy are not declared but correspond to wastewater treatment plants. Emission estimates referred to different European regions and from the SEF area are reported in Table B1.
If the inversion is run for each year of study period, an analogous source pattern is obtained. On the base of this result we extrapolated the E-PRTR figures to those years in which emission data were not available. The identification of the SEF region as the dominant source area in the study domain suggested a further analysis aimed at improving the emissions intensity estimate in the SEF area. An uncertainty of 100 % assumed for the a priori emissions used in previous test could indeed be insufficient to quantify emissions from a source much stronger with respect to the a priori emission field, such as the SEF region. Therefore, we repeated the inversion but with a dramatically increased uncertainty of 500 % in the SEF region. In the rest of the domain the uncertainty was maintained at 100 %. Hereinafter, we refer to this improved Bayesian inversion as BI-2. The emission map reported in Fig. B2 clearly confirms the SEF region as an emission hot spot responsible for the majority of the MCF emissions from the study domain.
The emission estimates for different European regions and from the SEF area as derived from the BI-2 analyses are reported in Table B1. For the sake of comparison, the emission estimates obtained using E-PRTR a priori field (reference inversion, RI) are reported as well. It is worth noting that the results obtained running the inversion with a homogeneous a priori emission field (BI-2) and with an a priori emission field based on the E-PRTR inventory (RI) differ by 18 %, suggesting high model robustness.