Top-down estimates of black carbon emissions at high latitudes using an atmospheric transport model and a Bayesian inversion framework

This paper presents the results of BC inversions at high northern latitudes (> 50 N) for the 2013–2015 period. A sensitivity analysis was performed to select the best representative species for BC and the best a priori emission dataset. The same model ensemble was used to assess the uncertainty of the a posteriori emissions of BC due to scavenging and removal and due to the use of different a priori emission inventory. A posteriori concentrations of BC simulated over Arctic regions were compared with independent observations from flight and ship campaigns showing, in all cases, smaller bias, which in turn witnesses the success of the inversion. The annual a posteriori emissions of BC at latitudes above 50 N were estimated as 560± 171 kt yr−1, significantly smaller than in ECLIPSEv5 (745 kt yr−1), which was used and the a priori information in the inversions of BC. The average relative uncertainty of the inversions was estimated to be 30 %. A posteriori emissions of BC in North America are driven by anthropogenic sources, while biomass burning appeared to be less significant as it is also confirmed by satellite products. In northern Europe, a posteriori emissions were estimated to be half compared to the a priori ones, with the highest releases to be in megacities and due to biomass burning in eastern Europe. The largest emissions of BC in Siberia were calculated along the transect between Yekaterinsburg and Chelyabinsk. The optimised emissions of BC were high close to the gas flaring regions in Russia and in western Canada (Alberta), where numerous power and oil and gas production industries operate. Flaring emissions in Nenets– Komi oblast (Russia) were estimated to be much lower than in the a priori emissions, while in Khanty-Mansiysk (Russia) they remained the same after the inversions of BC. Increased emissions at the borders between Russia and Mongolia are probably due to biomass burning in villages along the Trans-Siberian Railway. The maximum BC emissions in high northern latitudes (> 50 N) were calculated for summer months due to biomass burning and they are controlled by seasonal variations in Europe and Asia, while North America showed a much smaller variability.


Introduction
Light-absorbing species, such as black carbon (BC), are the main components of atmospheric particulate matter, affecting air quality, weather and climate.BC originates from the incomplete combustion of fossil fuels (primarily coal and diesel), from open high-temperature combustion of natural gas in the oil and gas fields (gas flaring), and from the burning of biomass and biofuels.BC particles affect cloud formation and precipitation as they act as cloud condensation nuclei in their hydrophilic form (Wang et al., 2016).BC is also a major driver of climate change, contributing to global warming with a radiative forcing at the top of the atmosphere ranging between 0.17 and 0.71 W m −2 (Bond et al., 2013;Myhre et al., 2013;Wang et al., 2014).BC deposited in Arctic snow surfaces in concentrations of up to 30 ng g −1 can reduce snow albedo by 1 %-3 % (Hegg et al., 2009) in fresh snow and up to 3 times more as snow ages and the BC particles become more concentrated (Clarke and Noone, 1985).Airborne BC also warms the air and reduces tropical cloudiness by absorbing the incoming solar radiation (Ackerman, 2000).It also reduces atmospheric visibility and increases N. Evangeliou et al.: Top-down estimates of black carbon emissions aerosol optical depth (Jinhuan and Liquan, 2000).From a health perspective, BC particles, generally being sub-micron in size, can penetrate into the lungs and cause pulmonary diseases (e.g.Wang et al., 2014).
To improve understanding about how BC affects climate and to develop effective policies to tackle BC's associated environmental problems requires accurate knowledge of the emissions and their spatio-temporal distribution.Most commonly, BC emission inventory datasets are built by "bottomup" approaches, which are based on activity data and emission factors and proxy information for spatial disaggregation, but these methods are considered to have large uncertainties (Cao et al., 2006).Numerous global or regional emission inventories of BC have been constructed previously (Bond et al., 2004;Schaap et al., 2004;Streets et al., 2003); nevertheless, emission uncertainties contribute significantly to the overall uncertainty of modelled concentrations of BC.Emission uncertainties affect regional or episodic simulations even more significantly, as in many cases emissions deviate from the annual mean.Such studies represent a useful tool to improve our understanding of the relationship between observed concentrations of BC and BC emissions.Furthermore, BC emissions have their most pronounced effect on the regional scale due to the relatively short atmospheric lifetime of BC (Hodnebrog et al., 2014;Samset et al., 2014),.
The relative differences between different emission inventories are largest for the high latitudes (AMAP, 2015), particularly in high-latitude Russia where emission information is poor.For this area, a new satellite-based high-resolution inventory showed that BC emissions from biomass burning (BB) might have been 3.5 times higher than emissions given in the Global Fire Emissions Database (GFEDv4) (Hao et al., 2016), if more realistic emission factors are used (May et al., 2014).Furthermore, new sources of BC in the same area have been identified recently.For example, emissions from gas flaring by the oil industry have been missing from most emission inventories and may be an important source of BC at high latitudes (Stohl et al., 2013).For instance, in 2008 Russia was responsible for nearly one-third of the gas flared globally (Elvidge et al., 2009).However, the gas flaring source is highly uncertain.For example, based on isotopic measurements, Winiger et al. (2017) reported recently that the contribution from gas flaring to BC measured at Tiksi in Siberia is lower than estimated by Stohl et al. (2013), while recently published bottom-up inventories (Huang et al., 2015;Huang and Fu, 2016) suggested even higher gas flaring emissions.Finally, Popovicheva et al. (2017) reported that one existing emission dataset of BC captured surface concentrations in the Russian Arctic quite efficiently.
In this study, we estimated the BC emissions at high northern latitudes using atmospheric observations of BC in a Bayesian inversion framework.Emissions were estimated for the region north of 50 • N because this is the region with the largest influence on Arctic surface concentrations (Klonecki, 2003;Stohl, 2006).We determine the emissions with monthly time resolution for the years 2013, 2014 and 2015.We first describe the observation data, the transport model and the Bayesian inversion technique used, as well as the prior emission information.We then assess the sensitivity of the transport model to different scavenging coefficients (below-cloud and in-cloud) for BC and to different emission inventories.We finally present optimised BC emissions, discuss these results in comparison with independent estimates and calculate the uncertainty of the inversions with respect to different scavenging parameters used for BC using four different prior emission datasets.

Observation network
Atmospheric observations of BC were retrieved from the World Data Centre for Aerosols (http://ebas.nilu.no,last access: 10 September 2016) and from the International Arctic Systems For Observing The Atmosphere (http://www.esrl.noaa.gov/psd/iasoa/, last access: 10 September 2016).An overview of the stations used in this paper can be found in Table 1 and Fig. 1a-c.The selected measurements were performed with different types of instruments that may differ substantially.When measurements are based on light absorption we refer to equivalent BC (EBC), while measurements based on thermal-optical methods refer to elemental carbon (EC) (Petzold et al., 2013).
At Alert (ALT), Appalachian (APP), Asprveten (ASP), Birkenes (BIR), East Trout Lake (ETL), South Great Planes (SGP), Steamboat Springs (COL), Trinidad Head (TRI) and Whistler (WHI) measurements were performed with particle soot absorption photometers (PSAPs).At Annaberg-Buchholz (ANB), Bösel (BOS), Cabauw Zijdeweg (CAB), Hyytiälä (HYY), Leipzig (LEI), Melpitz (MEL), Nepal Climate Observatory (NEP), Pallas (PAL), Ústí n.L.-mesto (ULM) and Waldhof (WLD) the particle light absorption coefficient was measured by multi-angle absorption photometers (MAAPs; Petzold and Schönlinner, 2004), which are in excellent agreement with other particle light absorption photometers such as a photoacoustic sensor (e.g.Müller et al., 2011).In the MAAP instrument, particles are continuously sampled on filter tape, with loaded spots subsequently analysed by Raman spectroscopy to derive the particle mass concentration of soot (Nordmann et al., 2013).The cut-off sizes of the different MAAP instruments varied between 1 and 10 µm.Continuous light absorption photometers (CLAPs, Model PSAP; 565 nm) were used at Utqiaġvik (formerly Barrow) (BAR), Bondville (BON), Gosan (GOS), Kpuszta (KPU), BEO Moussala (MOU) and Summit (SUM).Although these instruments were calibrated to measure the aerosol absorption coefficient, a previous study at this site revealed that a value of 10 m 2 g −1 is a reasonable conversion factor to determine the BC concentration (Gelencsér et al., 2000).Aethalometers were used at Tiksi (TIK) and Zeppelin (ZEP).All these stations measure the particle light absorption coefficients of different size fractions of the aerosol at wavelengths around 530-550 nm.Then the light absorption coefficients are converted to EBC mass concentrations under certain assumptions (Petzold et al., 2013).This is done externally for instruments such as MAAP, CLAP, PSAP etc. using a mass absorption efficiency of 10 m 2 g −1 (Bond and Bergstrom, 2006).For Aethalometers, the conversion is done internally by the instrument.All station measurements are routinely filtered to remove influence from local sources.

Source-receptor relationships (SRRs)
We used the Lagrangian Particle Dispersion Model (LPDM), FLEXible PARTicle dispersion model (FLEXPART) (Stohl et al., 1998(Stohl et al., , 2005) ) to model atmospheric transport.Using LPDMs to model particle or trace gas concentrations has several advantages over Eulerian models, namely that they can have quasi-infinite resolution and they are not subject to numerical diffusion.Thus they can provide better resolved source-receptor relationship (SRR) fields, which describe the relationship between the sensitivity of a "receptor" to a "source" element, as described by Seibert and Frank (2004).
SRRs for the lowest model level are often called footprint emission sensitivities or even just footprints.
SRRs were calculated using FLEXPART version 10 in a backwards mode (see Stohl et al., 2005) in which computational particles are released backward in time from the observation sites (receptors).When the number of observation sites is smaller than the number of unknown flux grid cells this mode is computationally more efficient than forward calculations.Furthermore, backward simulations can be initiated exactly at the measurement point without initial diffusion of information into a grid cell.This important advantage of LPDMs also facilitates high spatial resolution of the model output around the measurement sites.As meteorological input data, European Centre for Medium-Range Weather Forecasts operational meteorological analyses were used with 137 vertical levels and a horizontal resolution of 0.5 • × 0.5 • .Retroplumes were calculated at hourly intervals at each of the receptors.A total of 40 000 particles for each retroplume were released and followed 30 days backwards in time.This should be a sufficiently long time in order to include almost all contributions to BC concentration at the receptor given the atmospheric lifetime of BC (in the range of 2-10 days, Benkovitz et al., 2004;Koch and Hansen, 2005;Park et al., 2005;Textor et al., 2006).The treatment of scavenging is a major uncertainty for modelling BC (Browse et al., 2012).Therefore, an ensemble of 12 model simulations was performed, each with different BC tracers with different in-cloud and below-cloud scavenging properties (Table 2).This method allows the sensitivity of the SRRs (produced by FLEXPART) to scavenging to be quantified.Table 2 shows the different below-cloud and in-cloud scavenging parameters used within the model in the sensitivity runs.For all tracers, we assumed a logarithmic size distribution with an aerodynamic mean diameter of 0.25 µm, a logarithmic standard deviation of 0.3 and a particle density of 1500 kg m −3 (Long et al., 2013).The dry deposition scheme in FLEXPART is based on the resistance analogy (Slinn, 1982).The present version of the model uses the precipitation rate from ECMWF to determine below-cloud scavenging and the cloud liquid water and ice content, precipitation rate and cloud depth from ECMWF to calculate in-cloud scavenging (see Grythe et al., 2017).
The SRR at the lowest model layer (in seconds) (Fig. 1gi) can be multiplied with gridded emission fluxes from a BC emission inventory (in kg m −2 s −1 ) distributed over the layer depth (100 m).This gives the prior concentration of BC at the receptor point (in ng m −3 ).

Bayesian inversion
The inversion methodology used in the present study, FLEX-INVERT, is described fully in Thompson and Stohl (2014) and has been already used in studies of CH 4 , HFC-125, HFC-134a and SF 6 (Brunner et al., 2017;Thompson et al., 2015Thompson et al., , 2017)).Since atmospheric transport and deposition are lin-  (Giglio et al., 2013).(h) Monthly total (anthropogenic and biomass burning) BC emissions north of 50 • N from 2013 to 2015 from the four prior inventories used for the inversion.Coloured numbers correspond to total annual BC from each emission inventory.ear operations, they can be described as a Jacobian matrix of SRRs (H).The BC concentrations (y) can then be modelled given an estimate of the emissions (x) as follows: where ε is an error associated with model representation, such as the modelled transport and deposition as well as the measurements.Since H is generally not invertible (or may have no unique inverse), statistical optimisation methods are used, which require prior information for regularisation.According to Bayesian statistics, the problem can be expressed as the maximisation of the probability density function of the emissions given the prior information and observations and is equivalent to finding the minimum of the cost function: where B and R are the error covariance matrices for the prior emissions and the observations, respectively.The error in the observation space also accounts for model representation errors that are not related to the BC emissions.The emissions that minimise the cost function can be found by solving the first-order derivative of Eq. ( 2).Hence, the following equation can be derived for the most probable emissions, x (for details see, for example, Tarantola, 2005): In this study, the state vector contains the monthly unknown surface emissions on the grid of variable resolution (Fig. 1df) and has a resolution of between 1.0 • ×1.0 • and 8.0 • ×8.0 • .The total number of emission variables to be determined was 1422 for 2013, 1404 for 2014 and 1436 for 2015.The posterior error covariance matrix, A, is equivalent to the inverse of the second derivative of the cost function.However, to ac- A and B are rain and snow collection efficiencies for below-cloud scavenging; A i is the cloud condensation nuclei efficiency and B i the ice nuclei efficiency that are used in in-cloud scavenging following Grythe et al. (2017).2.0 1.0 0.45 0.10 BC9 0.2 0.2 0.90 0.90 BC10 1.0 1.0 0.90 0.20 BC11 2.0 1.0 0.45 0.45 BC12 1.0 1.0 0.45 0.00 count for the uncertainty in the scavenging parameters and different prior information, we instead conduct an ensemble of inversions to estimate the posterior uncertainty in order to account for the systematic errors.To do this, we conduct the inversion for BC represented by 12 different scavenging coefficients (see Table 2) and for four different prior emission datasets, and do this for each of the 3 years of our study (2013)(2014)(2015).The resulting model ensemble (12 × 4 = 48) for each year defines the posterior uncertainty due to scavenging and use of different a priori information (Sect.3.3).
Since negative values for the posterior emissions are mathematically possible but physically unlikely, we applied a subsequent inequality constraint on the emissions following the method of Thacker (2007).This is a truncated Gaussian approach in which inequality constraints are applied as errorfree "observations": where A is the posterior error covariance matrix, P is a matrix operator to select the variables that violate the inequality constraint, and c is a vector of the inequality constraint, which in this case is zero.The emissions were solved on an irregular grid, which has been optimized based on the SRRs to give higher resolution (1.0 • × 1.0 • ) in regions where there is strong contribution from emission sources to BC concentrations and lower (8.0 • × 8.0 • ) where there is a weak contribution (Stohl et al., 2009).Then, the results are interpolated onto a uniform grid of 1.0 • × 1.0 • resolution from 180 • W to 180 • E and 50 to 90 • N and are given at monthly time resolution for 2013, 2014 and 2015.To constrain emissions of BC monthly, a temporal correlation scale length between flux time steps equal to 90 days was set.

A priori emission information
In the present study, the emission inventories ECLIPSE (Evaluating the CLimate and Air Quality ImPacts of ShortlivEd Pollutants) version 5 (Klimont et al., 2017) (Wang et al., 2014) (available here: http://accent.aero.jussieu.fr/MACC_metadata.php, last access: 15 November 2016) were used as the prior emission estimates of BC (see Fig. 2).
The ECLIPSE emission inventory (Fig. 2a) accounts for waste burning, industrial combustion and processing, surface transportation that also includes power plants, energy conversion and extraction that also includes gas flaring, and residential and commercial combustion.
The HTAP_V2 dataset (Fig. 2b) consists of highresolution gridded emissions of BC based on nationally reported emissions combined with regional scientific inventories.It includes the sources of aviation, inland waterways and marine shipping, energy production other than electricity generation, industrial processes, solvent production and application, electricity generation, ground transport, buildings heating, cooling, equipment, and waste disposal or incineration.
The ACCMIP simulations use the BC emission inventory covering the historical period  provided by Lamarque et al. (2010), which is built for the climate model simulations in CMIP5 (Fig. 2c).Anthropogenic emissions are mainly based on Bond et al. (2004) but apply new emission factors.The year 2000 dataset was used for harmonisation with the future emissions determined by integrated assessment models (IAMs) for the four Representative Concentration Pathways (RCP4.5,RCP6, and RCP8.5).They include emissions from energy production and distribution, industry (combustion and non-combustion), transportation, maritime transport and aviation, residential and commercial combustion and solvent extraction, agricultural production, and waste treatment.
MACCity (Fig. 2d) was built as an extension of the historical emissions dataset of ACCMIP.It provides monthly averaged sectorial emissions for each year during the 1960-2010 period.This dataset was based on the decadal ACCMIP emissions for 1960-2000 and the 2005 and 2010 emissions provided by RCP 8.5.This scenario was chosen since it in- Figure 3. Taylor diagrams for the comparison of the prior (ECLIPSEv5) simulated concentrations with observations for all years (2013)(2014)(2015) for 12 BC species with different scavenging coefficients (Table 2).The radius indicates standard deviations normalised against the mean concentration (NSD); the azimuthal angle is the Pearson correlation coefficient, while the normalised (against observation) root mean square error (nRMSE) in the simulated concentrations is proportional to the distance from the point on the x axis identified as "reference" (grey contours).
cluded some information on recent emissions at the regional scale in Europe and North America.The emission sectors are consistent with Lamarque et al. (2010).
Emissions from biomass burning were adopted from the Global Fire Emissions Database, Version 4 (GFEDv4) (Giglio et al., 2013), and implemented to each of the four emission inventories for 2013, 2014 and 2015.Emissions from gas flaring are only included in ECLIPSEv5 inventory.

Sensitivity to scavenging and selection of the best representative species for BC
The comparison of the simulated and observed concentrations for the 12 different BC tracers (see Table 2) is shown in Taylor diagrams in Fig. 3   Correlations of modelled and observed surface concentrations of BC were high (> 0.5) only at stations ALT, MEL and LEI, which present low normalised standard deviation (NSD) values (< 1) and low normalised root mean square error (nRMSE) values.All NSD values were below 1.5 except at TIK and ULM stations (see Figs. 3 and S1).In general, dispersion models fail to reproduce BC concentrations close to TIK station (Eckhardt et al., 2015;Evangeliou et al., 2016), as the station has been reported to receive pollution from local anthropogenic sources (Asmi et al., 2016).ULM station is located on the border between Germany and the Czech Republic and was shown previously to be strongly affected by BC emissions from residential combustion sources (Schladitz et al., 2015).The model-observation mismatches ([model − observations] / observations) due to the use of 12 different species for BC can be seen in Fig. S2 for the years 2013, 2014 and 2015.These values are average concentrations from the use of the four different emissions inventories (ECLIPSEv5, ACCMIPv5, EDGAR_HTAPv2.2and MAC-City).The extreme perturbation of the scavenging coefficients of BC caused an average relative model-observation mismatch (normalised against observations) of about 39 % in 2013 (Fig. S2) at all stations.
Similar to 2013, 12 species with different scavenging parameters were used for BC following Table 2 in 2014 and the comparison with observations is shown in Taylor diagrams in Fig. 3 using ECLIPSEv5 emissions and in Fig. S1 using AC-CMIPv5, EDGAR_ HTAPv2.2 and MACCity for the common stations.The comparison of surface simulated concentrations with observations showed NSD values above 1, high nRMSE values and correlation coefficients below 0.5 in at most of the stations.The main difference from year 2013 is that the model-observation mismatches for the surface concentrations of the 12 BC species was estimated to be 32 % in 2014 (Fig. S2), in contrast to 39 % in 2013.The same deficiency of the model to capture the spring and summer concentrations of BC was observed.The calculated mismatches were very low in at most of the lower latitude stations and increased towards the remote Arctic ones (Fig. S2).
Finally, in 2015 the comparison of surface concentrations for each of the 12 different BC species using the four different datasets (ECLIPSEv5, ACCMIPv5, EDGAR_HTAPv2.2and MACCity) with observations showed again the same pattern as in the previous years, with most of the NSD values being above unity, high nRMSE values and low Pearson coeffi- cients (Figs. 3 and S1).The model-observation mismatches of BC concentrations (Fig. S2) were estimated to be as high as 43 % for the stations where full measurements existed for the 3 years of the study (2013)(2014)(2015).Like in the previous years, the model failed to reproduce surface concentrations of BC at some of the remote stations of the Arctic.
We used the NMSE (normalised mean square error) to select the most representative BC tracer species.The NMSE is an estimator of the overall deviations between predicted and measured values.It is defined as follows: where O i and P i are observed and predicted concentrations, N is the number of observations for which we assess the predicted values, and the overbar indicates the mean over the number of observations for O i and P i .Contrary to the relative mismatches, in the NMSE the squared deviations (absolute values) are summed instead of the differences.For this reason, the NMSE generally shows the most striking differences among models.NMSE is a highly selective statistical quantity that can give large differences between models that perform similarly for other statistical measures.The lower the NMSE value, the better the performance of the model.On the other hand, high NMSE values do not necessarily mean that a model is completely wrong as the errors could be due to shifts in time and/or space.Moreover, it must be pointed out that NMSE is sensitive to outliers (Poli and Cirillo, 1993).
The calculated monthly average NMSE values for the 12 species using ECLIPSEv5 as the emission input can be seen in Fig. 4 for the years 2013-2015.The different scavenging coefficients used did not create a large variability in the monthly BC concentrations.This was caused due to the small perturbation of the scavenging coefficients.A more drastic change of wet scavenging would have caused BC concentrations to change more, as wet scavenging dominates removal and deposition of BC by approximately 80 %.However, considering that the selection of the best representative species for BC was a top priority, an effort to set realistic values to the scavenging coefficients was made.The best performance for the majority of the stations examined and most months was obtained for species 1, 2 and 10 (see Table 2).In terms of model response over the Arctic stations, a better performance was achieved for species 1 than for the other two.Therefore, we have chosen species 1 as our reference species for all subsequent analyses and the inversions.It should be noted here that the same test was performed using ACCMIPv5, EDGAR_HTAPv2.2and MACCity emissions.Although the results were worse, the best-performing species for BC were again 1, 2 and 10.

Sensitivity to different prior information and selection of the best prior emission inventory
In this section we assess the impact of using the different prior emission inventories for BC and select the most ap- propriate one for our BC inversions.For this analysis, the best-performing species 1 for BC (see Table 2) was chosen and the monthly relative model-observation mismatches (([model − observations]/observations) for all stations and years separated were calculated using all four inventories and are depicted in Fig. S3.
The largest monthly relative model-observation mismatches for the a priori simulated concentrations of BC in 2013 were calculated for stations located close to 50 • N (BOS, CAB MEL, LEI, ULM, ANB).The average modelobservation mismatch for all stations was 15 % for 2013.Similar results were found for 2014 for the prior simulated BC concentrations with the largest relative mismatches recorded at mid-latitude stations where the BC concentrations were very high due to large anthropogenic emissions (BOS, CAB, MEL, LEI).On average, the relative modelobservation mismatch was as high as 23 % for the year 2014.Finally, in 2015, the highest monthly relative mismatches of the a priori BC concentrations were again estimated for the stations of high anthropogenic influence (CAB, MEL, LEI, WLD).The average relative model-observation mismatch in 2015 was only 19 %, much lower than all previous years.The fact that all prior emission datasets used failed to reproduce the observations in central Europe during all years studied (2013, 2014 and 2015), whereas other stations at mid-latitudes were reproduced well, might imply either missing sources or highly uncertain measurements (Fig. S3).The use of different emission dataset changes simulated concentrations by a maximum of 23 %.
Normalised mean square error (NMSE) values calculated for each of the four emission inventories were very low at the majority of the stations for which data existed in all the years of study (ZEP, SUM, TIK, BAR, MEL and LEI), when ECLIPSEv5 emissions were used.In contrast, at PAL all emission datasets performed well (Fig. 5).The observations of BC concentrations at Arctic stations were better reproduced in simulations using the ECLIPSEv5 than with any other inventories examined.Law and Stohl (2007) have documented that these elevated BC concentrations are caused by anthropogenic emissions.Black carbon concentrations at TIK are not well simulated for reasons given in Sect.3.1.

Optimised (a posteriori) emissions of BC and associated uncertainty
The optimised annual emissions of BC together with the associated posterior gridded uncertainty and the difference between posterior and prior emissions averaged for the 2013-2015 period can be seen in Fig. 6.The posterior emissions are presented for the best-performing species (species 1) of BC and the best prior emissions inventory (ECLIPSEv5).The to- The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table 2. Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets (Sect.3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.
tal posterior uncertainty was calculated as the standard deviation of the posterior emissions calculated for the 12 BC species with different scavenging coefficients for four different emission datasets as prior information for each of the 3 years (12 × 4 = 48 inversions; see Sect.2.3).The total uncertainty is a propagation of the deposition uncertainty (represented by the posterior emissions using 12 perturbed BC species with different scavenging coefficients) and the uncertainty due to the use of different prior information (represented by the posterior emissions using the four different emission datasets).The optimised emissions show some constant hotspot areas that persist throughout all 3 years and which are attributed to anthropogenic BC emissions.For instance, emissions in the Nenets-Komi region close to the Yamal peninsula in Russia or in Khanty Mansiysk region of northwestern Siberia have been reported to originate to a large extent from gas flaring (Popovicheva et al., 2017;Stohl et al., 2013;Winiger et al., 2017).Other areas that are characterised by large anthropogenic emissions are in western Canada (Alberta), where more than 100 power industries burn fossil fuels and more than 50 oil and gas production and oil-refining industries operate.In addition, one of the largest oil sands deposits is found in northern Alberta and in the McMurray area, which contains about 168 billion barrels of oil (Heins, 2000).Cheng et al. (2018) found high concentrations of BC (more than 1000 ng m −3 ) in the Canadian oil sands region at altitudes of up to 2 km during a flight campaign.
The optimised BC emissions in North America for the 2013-2015 period were between 149 and 193 kt yr −1 (average ±1σ error: 174 ± 58 kt yr −1 ), in the same order as the prior emissions in ECLIPSEv5 (148-182 kt yr −1 ) and slightly higher than ACCMIPv5, EDGAR_HTAPv2.2and MACCity (116-150 kt yr −1 ).In northern Europe we estimated that 124-238 kt yr −1 of BC was released (average ± 1σ error: 170 ± 59 kt yr −1 ), which is less than half the ECLIPSEv5 emissions (352-381 kt yr −1 ), about 35 % lower than the ACCMIPv5 and MACCity emissions (241-256 kt yr −1 ), and in the same order as the EDGAR_HTAPv2.2emissions (163-175 kt yr −1 ).Posterior BC emissions were higher in northern Siberia for the 3year period (130-291 kt yr −1 , average ±1σ error: 217 ± 69 kt yr −1 ) compared with ECLIPSEv5 ( 187 Shevchenko et al. (2016).The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table 2. Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets (Sect. 3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.
ACCMIPv5 (127-178 kt yr −1 ), EDGAR_HTAPv2.2(108-159 kt yr −1 ) or MACCity (129-179 kt yr −1 ).Larger changes in the BC emissions were calculated in Russian territories that are known to be important gas flaring sources (Stohl et al., 2013).BC emissions in the Nenets-Komi oblast were between 14 and 17 kt yr −1 (average ±1σ error: 15 ± 5 kt yr −1 ), about 40 % lower than the respective emissions in ECLIP-SEv5 (≈ 25 kt yr −1 ), the only prior dataset that took gas flaring into account there.This could be due to the decreasing magnitude of the flaring emissions in the last few years (see Huang and Fu, 2016).Finally, in Khanty-Mansiysk BC emissions were 28-37 kt yr −1 (average ±1σ error: 32 ± 8 kt yr −1 ) compared to 25 kt yr −1 in ECLIPSEv5, whereas in the other datasets that do not include BC emissions due to flaring, BC emissions were negligible.However, the posterior Khanty-Mansiysk emissions are shifted further east compared to the prior.
The annual posterior BC emissions at latitudes above 50 • N were estimated as 560 ± 171 kt yr −1 averaged for the 2013-2015 time period (523 ± 92 kt yr −1 in 2013, 608 ± 104 kt yr −1 in 2014 and 549 ± 100 kt yr −1 in 2015, respectively).For the same area and period, BC emissions in ECLIPSEv5 were 745 kt yr −1 , in ACCMIPv5 533 kt yr −1 and in EDGAR_HTAPv2.2437 kt yr −1 , while in MACCity they were 538 kt yr −1 .The annual posterior absolute uncertainty can be seen in Fig. 6b.As it was explained before, this uncertainty is a combination of the uncertainty due to scavenging and due to the use of different prior information in the inversions of BC.Averaged over the period 2013-2015, the relative uncertainty of the inversion was estimated to be 30 %.The uncertainty due to different scavenging coefficients in the BC species used was 25 %, while the uncer-tainty due to the use of different prior emissions was only 5 %.

Validation of the posterior emissions of BC
The concentrations of BC at eight measurement stations simulated with the posterior (optimised) BC emissions can be seen in Fig. S4.As expected, BC concentrations match the observations significantly better than using any of the a priori datasets with correlation coefficients above 0.6 for most of the stations.At the same time, NSD values were close to unity or lower and the nRMSE values were below 1.5 at most of the stations shown in Fig. S4.However, the comparison to observations included in the inversion is not a sufficient indicator of the inversion's performance, as the inversion is designed to reduce the model-observation mismatches.The magnitude of the posterior reduction of the model mismatch to the observations is partly determined by the weighting given to the observations relative to the prior emissions.A much better performance indicator is the comparison of the posterior concentrations with observations that were not included in the inversion (independent observations).
For this reason, we compared posterior BC concentrations with observations from the ACCACIA (Aerosol-Cloud Coupling and Climate Interactions in the Arctic) flight campaign, which was conducted near Zeppelin station, Ny-Ålesund, for 3 days in March 2013 (Sinha et al., 2017).This campaign was chosen because it was conducted during 1 year for which inversion results are available (2013).The results are shown in Fig. 7 for the prior simulated concentrations of BC using four different emission datasets (ACCMIPv5, ECLIPSEv5, EDGAR_HTAPv2.2and MACCity) and the posterior simu-  2015) with observations from a ship campaign in the Russian Arctic in 2015 adopted from Popovicheva et al. (2017).The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table 2. Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets (Sect.3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.lated BC concentrations.In all profiles, use of the optimised BC emissions results in a better agreement between modelled concentrations and observations compared to the prior simulated BC concentrations, while the RMSE (not normalised) values decrease substantially.However, the Pearson's correlation coefficients were below 0.5.
To assess the performance of the inversions of BC in 2014, we used an independent dataset from a ship campaign that took place in the North Atlantic and Baltic Seas in June and August 2014 (Fig. 8), provided by Shevchenko et al. (2016).Although the measurements may sometimes be affected by the ship's exhaust, the posterior RMSE was 34 % lower than the average RMSE using four different a priori emission datasets (ACCMIPv5, ECLIPSEv5, EDGAR_HTAPv2.2and MACCity), supporting the view that the inversion improved the emissions for 2014.
To validate the 2015 inversions of BC, measurements from a ship campaign over the Russian Arctic were used (Popovicheva et al., 2017) and the results are shown in Fig. 9.The cruise started from the port of Arkhangelsk in the northwestern European Russia, reached the Bolshevik Island in the higher Russian Arctic and returned following more or less the same pathway.The calculated RMSE of the posterior BC concentrations with the measurements taken during the cruise was about 10 % lower that the respective RMSE from the prior simulated concentrations of BC (average for all prior simulated emissions).This shows that the optimised emissions improved BC concentrations over the Russian Arctic.Some episodic peaks of BC throughout the ship cruise, however, were poorly captured.

BC emissions in North America
The spatial distribution of the optimised BC emissions in North America averaged for the 3-year period is depicted in Fig. 10 and the annual posterior emissions for 2013, 2014 and 2015 are shown in Fig. S5. Figure 10b shows the differences between posterior and prior emissions (ECLIPSEv5) and highlights the biggest emission changes compared to the a priori dataset.
The most characteristic locations of sources between 2013 and 2015 lie in Alberta, where most of the large oilproducing industries operate (Figs. 10 and S5).The highest emission source was located in 60 • N, 135 • W in 2013 and 2015, but not in 2014.This spot corresponds to the location of Whitehorse, which is the capital and only city of Yukon, and the largest city in northern Canada.The area contains mining activities (mainly for gold) and three natural gas wells, while biomass in the form of cordwood and pellets is used for space heating (Yukon Government, 2018).The fact that near-zero BC emissions were calculated in Whitehorse in 2014 might be due to the lack of available measurements in North America, which in turn results in poorly constrained posterior BC emissions.Another similar hotspot area that is more intense in 2013 and 2015, but not in 2014, is located in Yellowknife, north of Great Slave Lake (62.5 • N, 115 • W, Fig. S5).The city is known for gold and diamond mining and an oil-driven power plant (Northwest Territories Power Corporation, https://www.ntpc.com, last access: 5 March 2018).Finally, another characteristic hotspot emission region of BC is seen southeast of Lake Athabasca (57 • N, 108 • W, Fig. 10).Uranium mines are located in this region.These mines use diesel generators, diesel trucks and other diesel-powered machinery.Exactly in this location, the Visible Infrared Imaging Radiometer Suite (VIIRS) showed relatively strong night-time light sources (https://www.lightpollutionmap.info/#zoom=5&lat=8255540&lon=-11864816&layers=B0FFFTFFFT, last access: 5 March 2018).

BC emissions in northern Europe
The posterior BC emissions in northern Europe averaged for the period 2013-2015 can be seen in Fig. 11 together with the difference between prior (ECLIPSEv5) and posterior BC emissions, while the posterior emissions for each individual year are shown in Fig. S6.The location of the gas flaring facilities are also presented in the same figures together with vegetation fires from the FEINE (Fire Emission Inventory -northern Eurasia) inventory (Hao et al., 2016).The latter combines the MODIS thermal anomaly products (MOD14 and MYD14) and the MODIS top-of-theatmosphere-calibrated reflectance product (MOD02) to map and date burn scars that are screened for false detections.29 February 2017), grey points show the power industries that operate using fossil fuels and oil and gas production and oil-refining industries adopted from Industry About (https://www.industryabout.com/canada-industrial-map,last access: 29 ebruary 2017), and dark green points show active fires from MODIS.
Land cover classification of burned areas are taken from the MODIS land cover change product (MOD12) (Friedl et al., 2010).This dataset is considered to be more realistic than GFED4 due to the emission factors used for BC (May et al., 2014) and the different approach to burned area calculation (see Hao et al., 2016).
The highest posterior BC emissions are calculated for the Moscow megacity at 55 • N, 37.5 • E, Berlin at 52 • N, 14 • E, Warsaw at 52 • N, 21 • E, Kyiv at 50 • N, 30 • E and Saint Petersburg at 60 • N, 30 • E, while London is slightly misplaced to the west (Fig. 11).The Scandinavian countries have the lowest emissions, although domestic heating there can also  be important (Andersen and Jespersen, 2016).The difference between prior and posterior emissions show that vegetation fires have a large impact on the BC emissions, especially in eastern Europe.In particular in 2015, the inversion produces large emission increases exactly where a large number of fire hot spots were found (see Fig. S6).

BC emissions in northern Siberia
Figure 12 illustrates the average posterior BC emissions in western Siberia for the 2013-2015 period together with the difference between the prior (ECLIPSEv5) and the posterior BC emissions, while Fig. S7 shows the respective BC emis-sions for each year individually, together with the flaring facilities and the vegetation fires, similarly to the previous section.
The prior BC emissions from flaring in Nenets-Komi oblast are confirmed by the inversion, although the emissions are shifted further east, while the flaring emissions in Khanty-Mansiysk are probably underestimated in ECLIP-SEv5 (see also Table 3).Krai, Russia, due to the Norilsk Mining and Smelting Factory extracting coal and ores.Furthermore, an increase in posterior BC emissions was estimated across the line that connects some important Russian cities (Yekaterinsburg to Chelyabinsk, 55 • N, 60 • E).These cities have been reported to contribute large amounts of BC mainly from transportation (see Evangeliou et al., 2018).Another hotspot exists at 58  12).
In the western part of Siberia, there are numerous sources of average or low intensity.However, there no known anthropogenic sources there.At the lowest part of the inversion domain, in the borders of Russia with Mongolia, the posterior emissions showed a large increase (Fig. 6).These emissions are prevalent along the Trans-Siberian Railway.Human activities in the villages along the railway have been highlighted to be the major cause of the fires there.

Seasonal variability of BC emissions
The monthly optimised BC emissions are shown in Fig. 13 for the 3 years of study (2013)(2014)(2015) for the entire area north of 50 • N, and separately for areas north of 50 • N in North America, northern Europe, northern Siberia, Nenets-Komi oblast (Russia) and Khanty-Mansiysk (Russia).The last two regions are known to have large emissions from gas flaring.In the same figure the prior emissions from ECLIPSEv5, EDGAR_HTAPv2.2,ACCMIPv5 and MACCity are plotted for comparison.
The total posterior BC emissions (> 50 • N) show large seasonal variation (Fig. 13a).The maximum emissions were calculated for summer months (July in 2013 and June in 2014 and 2015).In these months large emissions from biomass burning have been reported both in GFED4 (see burned area in Giglio et al., 2013), and in FEINE (Hao et al., 2016).Separating the inversion domain into continental regions reveals where biomass burning is important.For instance, in North America (Fig. 13b) our optimised emissions show a significantly smaller variability, although GFED4 is included in all the prior emission datasets and shows a large emission peak for BC in summer, implying that fires are important.This was not the case for northern Europe (Fig. 13c), where the largest seasonal BC emissions were found in July for 2013 and in May for 2015, while in 2014 the largest peak appeared in April.This is not seen in the prior emission datasets, which show weak monthly variation.The largest seasonal variations were calculated for northern Siberia (Fig. 13d) and BC emissions there control the overall seasonal pattern for the total optimised BC emissions (> 50 • N).A large month-tomonth variability was estimated in the Nenets-Komi oblast (Fig. 13e), but no clear seasonal pattern was found.Finally, the largest monthly BC emissions in Khanty-Mansiysk oblast of Russia (Fig. 13f) were calculated in April for 2013, July for 2014 and June for 2015, showing that a large share of the BC emissions in this region originates from biomass burning since the region is located at mid-latitudes (60 • N-65 • N) and is vulnerable to open fires.

Conclusions
We have optimised BC emissions at high northern latitudes (> 50 • N) for the 2013-2015 period using a Bayesian inversion tool, an atmospheric transport model and network of continuous measurements of BC.We performed a sensitivity study to assess the best representative species for BC according to the efficiency of in-cloud and below-cloud scavenging, and the best representative emission inventory to be used as the prior information for our inversion.
The perturbation of scavenging coefficients for BC in the simulated concentrations creates a relative modelobservation mismatch of 32 %-43 % for the 3 years of study, whereas the use of different emission inventories has a less significant effect on the simulated concentrations, showing a relative model-observation mismatch of 15 %-23 %.
The posterior BC emissions show characteristic hotspots throughout all 3 years in the Nenets-Komi region close to the Yamal peninsula in Russia or in Khanty Mansiysk region of northwestern Siberia, where gas flaring facilities are located, and in western Canada (Alberta), where more than 150 power and oil-gas production industries operate.The annual posterior BC emissions at latitudes above 50 • N were estimated as 560 ± 171 kt yr −1 , significantly smaller than in ECLIPSEv5 (745 kt yr −1 ), which was used in the prior information in the inversions of BC (best representative emission dataset).
The uncertainty of the inversions was assessed using a model ensemble represented by 12 different scavenging coefficients for BC and four different prior emission datasets (12 × 4 = 48) for each of the 3 years of our study.We calculate a relative uncertainty of the inversion of 30 % for the 3 years of our study.
The posterior simulated concentrations of BC showed a better agreement with independent observations adopted from flight and ship campaigns over the Arctic presenting, in all cases, up to three times lower RMSE values.
In North America, the posterior emissions were found to be similar to the a priori ones driven by anthropogenic sources, while biomass burning appeared to be insignificant.This was confirmed by satellite products that showed weak existence of active fire hotspots.
In northern Europe, posterior emissions were estimated to be half compared to the prior ones, with the highest releases to be in megacities and due to biomass burning in eastern Europe.
Finally, in northern Siberia the larger emissions were calculated along the transect between Yekaterinsburg and Chelyabinsk, while flaring in Nenets-Komi oblast is probably overestimated in the a priori emissions.Increased emissions in the borders between Russia and Mongolia are probably due to biomass burning in villages along the Trans-Siberian Railway.
Data availability.All data generated for the present publication are stored on NIRD (https://www.uio.no/english/services/it/research/storage/nird-sigma.html,last access: 19 October 2018) (project NS9419K) and can be obtained from the corresponding author upon request.
The Supplement related to this article is available online at: https://doi.org/10.5194/acp-18-15307-2018supplementAuthor contributions.NE performed the simulations and analyses and wrote the paper.RLT helped in the adaptation of FLEXINVERT for BC and commented on the paper.SE helped in the implementation of the experiments, and AS coordinated, commented on and wrote parts of the paper.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Observation network used for the present inversion (a, b, c), and variable-resolution grid used for the inversion (d, e, f) also showing the location of the observation sites (red stars) for the 2013-2015 period.Sensitivity to the surface emissions (i.e. the footprint emission sensitivity or equivalent source-receptor relationship) integrated over all observation sites and all time steps (g, h, i) for the years 2013, 2014 and 2015 (units of log(ns)).

Figure 2 .
Figure 2. (a-d) Anthropogenic emissions of BC in the inversion domain (> 50 • N) from ECLIPSEv5, EDGAR_ HTAPv2.2,ACCMIPv5 and MACCity (anthropogenic emissions are assumed to be constant throughout every year).(e-g) Biomass burning emissions from GFED4 for 2013, 2014 and 2015(Giglio et al., 2013).(h) Monthly total (anthropogenic and biomass burning) BC emissions north of 50 • N from 2013 to 2015 from the four prior inventories used for the inversion.Coloured numbers correspond to total annual BC from each emission inventory.

Figure 4 .
Figure 4. Monthly average NMSE values due to use of 12 different BC species defined in Table2for the eight stations with complete data in the period 2013-2015.The annual mean is denoted as "ave").

Figure 5 .
Figure 5. Monthly average NMSE values due to use of different emission inventories (ECLIPSEv5, ACCMIPv5, EDGAR_HTAPv2.2,MACCity) for the eight stations with complete data in the period 2013-2015.The annual mean is denoted as "ave").

Figure 6 .
Figure 6.(a) Annual posterior emissions of BC in areas > 50 • N averaged for the period 2013-2015.(b) Average posterior uncertainty due to scavenging and use of different prior emissions for the same period.(c) Difference between posterior and prior emissions for 2013-2015.

Figure 7 .
Figure 7.Comparison of prior (ECLIPSEv5, ACCMIPv5, EDGAR HTAPv2.2 and MACCity) and posterior simulated concentrations of BC with observations from the ACCACIA flight campaign near Zeppelin station, Ny-Ålesund, in 2013, adopted from Sinha et al. (2017).The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table2.Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets(Sect.3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.

Figure 8 .
Figure8.Comparison of prior (ECLIPSEv5, ACCMIPv5, EDGAR HTAPv2.2 and MACCity) and posterior simulated concentrations of BC with observations from a ship campaign in North Atlantic and Baltic Seas in 2014 adopted fromShevchenko et al. (2016).The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table2.Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets(Sect.3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.

Figure 9 .
Figure 9.Comparison of prior (ECLIPSEv5, ACCMIPv5, EDGAR HTAPv2.2 and MACCity) and posterior simulated concentrations of BC (2015) with observations from a ship campaign in the Russian Arctic in 2015 adopted fromPopovicheva et al. (2017).The variability of the prior concentrations (shaded area) was calculated as the standard deviation of BC concentrations from the 12 species with different scavenging coefficients as shown in Table2.Uncertainties of the posterior concentrations are due to scavenging and use of four different a priori datasets(Sect.3.4).RMSE values are computed for ECLIPSEv5 concentrations, all prior concentrations (average) and posterior simulated BC concentrations.

Figure 10 .
Figure 10.(a) Optimised emissions of BC in North America (western Canada) averaged over the 2013-2015 period.(b)Difference between a posteriori and a priori emissions of BC (ECLIPSEv5 was used as the prior).Magenta points on the map denote the gas flaring industries from the Global Gas Flaring Reduction Partnership (GGFR) (http://www.worldbank.org/en/programs/gasflaringreduction,last access: 29 February 2017), grey points show the power industries that operate using fossil fuels and oil and gas production and oil-refining industries adopted from Industry About (https://www.industryabout.com/canada-industrial-map,last access: 29 ebruary 2017), and dark green points show active fires from MODIS.

Figure 11 .
Figure 11.(a) Optimised emissions of BC in northern Europe averaged over the 2013-2015 period.(b) Difference between a posteriori and a priori emissions of BC (ECLIPSEv5 was used as the prior).Magenta points on the map indicate the gas flaring industries from the Global Gas Flaring Reduction Partnership (GGFR) (http://www.worldbank.org/en/programs/gasflaringreduction,last access: 29 February 2017), while dark green points show the vegetation fires adopted from Hao et al. (2016).

Figure 12 .
Figure 12.(a) Optimised emissions of BC in western Siberia averaged for the 2013-2015 period.(b) Difference between a posteriori and a priori emissions of BC (ECLIPSEv5 was used as the prior).Magenta points on the map indicate the gas flaring industries from the Global Gas Flaring Reduction Partnership (GGFR) (http://www.worldbank.org/en/programs/gasflaringreduction,last access: 29 February 2017), while dark green points show the vegetation fires adopted from Hao et al. (2016).
Vegetation fires are shown to correlate well with BC emissions for 2013 (60-70 • N) and 2014 (50-60 • N) (Fig. S7), but not in 2015.Hotspots of high emissions were found in Dudinka (70 • N, 88 • E), a town on the Yenisei River and the administrative centre of Taymyrsky Dolgano-Nenetsky District of Krasnoyarsk

Figure 13 .
Figure 13.Monthly posterior emissions of BC shown for all regions located (a) > 50 • N, (b) in North America (> 50 • N), (c) in northern Europe (> 50 • N), (d) in northern Siberia (> 50 • N), (e) in Nenets-Komi oblast (> 50 • N, Russia) and (f) in Khanty-Mansiysk oblast (> 50 • N, Russia) for the 2013-2015 period.Monthly prior emissions of BC from ECLIPSEv5, EDGAR_HTAPv2.2,ACCMIPv5 and MACCity emissions inventories are also shown for the same regions and time period.The uncertainty of the posterior emissions of BC stems from the use of different scavenging coefficients and different prior emission datasets (see Sect. 3.3).

Table 1 .
Observation sites used for the inversions (the altitude indicates the sampling height in metres above sea level).

Table 2 .
Different scavenging parameters of below-cloud and incloud scavenging used in the ensemble model simulations for BC.