Source apportionment of black carbon and combustion-related CO₂ for the determination of source-specific emission factors

Black carbon (BC) aerosol typically has two major sources in the urban environment: traffic and domestic biomass burning, which has a significant contribution to urban air pollution during the heating season. Traffic emissions have been widely studied by both laboratory experiments (individual vehicle emission) and real-world measurement campaigns (fleet emission). However, emission information from biomass burning is limited, especially an insufficiency of experimental results from real-world studies. In this work, the black carbon burden in the urban atmosphere was apportioned to fossil fuel (FF) and biomass burning (BB) related components using the Aethalometer source apportionment model. Applying the BC source apportionment information, the combustion-related CO₂ was apportioned by multilinear regression analysis, supposing that both CO₂ components should be correlated with their corresponding BC component. The combination of the Aethalometer model with the multilinear regression analysis (AM-MLR) provided the source-specific emission ratios (ERs) as the slopes of the corresponding BC–CO₂ regressions. Based on the ER values, the source-specific emission factors (EFs) were determined using the carbon content of the corresponding fuel. The analysis was carried out on a 3-month-long BC and CO₂ dataset collected at three monitoring locations in Ljubljana, Slovenia, between December 2019 and March 2020. The measured mean site-specific concentration values were in the 3560–4830 ng m⁻³ and 458–472 ppm ranges for BC and CO₂, respectively. The determined average EFs for BC were 0.39 and 0.16 g(kg fuel)⁻¹ for traffic and biomass burning, respectively. It was also concluded that the traffic-related BC component dominates the black carbon concentration (55 %–64 % depending on the location), while heating has the major share in the combustion-related CO₂ (53 %–62 % depending on the location). The method gave essential information on the source-specific emission factors of BC and CO₂, enabling better characterization of urban anthropogenic emissions and the respective measures that may change the anthropogenic emission fingerprint.

1 Introduction

Biomass burning (BB) is a significant source of black carbon (BC), brown carbon (BrC), and organic particulate matter, creating a contribution to climate change (Myhre et al., 2013; Tomlin, 2021) and a severe risk to human health (Naeher et al., 2007; Janssen et al., 2011; Sigsgaard et al., 2015; Chen et al., 2017; Brown et al., 2020; Karanasiou et al., 2021). Global hotspots of BB are associated with extensive and persistent wildfires (e.g., deliberate forest burning in Amazonia and Indonesia, accidental forest and savanna fires in central Africa, North America, the Mediterranean basin, and Siberia) (see, e.g., Val Martin et al., 2006; Giglio et al., 2013; Smirnov et al., 2015; Chiloane et al., 2017; Healy et al., 2019; Reddington et al., 2019). On the other hand, emission from domestic wood combustion for the purpose of space heating, water boiling, or cooking significantly contributes to the BB emission as well, especially in locations of high population density and reduced ventilation (Karagulian et al., 2015; Klimont et al., 2017; Mitchell et al., 2017).

Wood combustion is an important energy source, even in well-developed countries, where its emissions add to traffic-related air pollution. The share of wood combustion in the total European energy budget is expected to increase dramatically due to the current energy crisis that most European countries face nowadays.

The emission characteristics of BB differ from those of internal combustion of fossil fuel (FF), where the combustion is more complete. Consequently, FF combustion emits less CO, particulate matter, and organic compounds per unit of fuel mass while having higher NO_X emissions compared to BB due to the higher combustion temperature and excess of air (see EEA, 2019: 1.A.3.b. versus 1.A.4.a–b.).

Black carbon is a dominant form of particulate matter emitted from fossil fuel combustion. Diesel engines (before the Euro 5 legislation standard) emit more than 80% of the particle mass (PM) as BC (EEA, 2019: 1.A.3.b.). Since diesel vehicles dominate the European vehicle fleet (Cooper, 2020), the high traffic-related BC emission poses significant air quality problems in cities, which is complemented by the BB emission during the heating season.

Due to their harmful health effects, BC emissions of diesel engines have been studied intensively worldwide for a long time (see the original work of Hansen and Rosen, 1990). The BC emission factors have been determined by numerous studies based on laboratory chassis dynamometer tests (Alves et al., 2015; Park et al., 2020) or real-world on-road measurements using either the chasing method (Wang et al., 2012; Ježek et al., 2015; Zavala et al., 2017) or on-board tailpipe measurements by the PEMS (Portable Emission Measurement System) (Zheng et al., 2015; Giechaskiel et al., 2019). These tests refer to the emission factors (EFs, emitted pollutant per kilogram of fuel or kilometer) of individual vehicles and do not reflect the emission of the entire vehicle fleet. By contrast, roadside monitoring offers the opportunity to measure a statistically significant number of vehicles. These measurements are usually carried out in tunnels (Sánches-Ccoyllo, 2005; Ban-Weiss et al., 2009; Brimblecombe et al., 2015; Blanco-Alegre et al., 2020), where elevated pollution concentration levels and negligible interference of other combustion sources (like wood burning) can be ensured. In these studies, the EF calculation is usually based on the carbon-balance method (Brimblecombe et al., 2015), when the plume CO₂ increment is used to determine the burnt fuel mass.

In contrast, emissions from biomass burning are not controlled nearly as strictly as from mobile sources. Some studies have investigated specific combustion appliances, providing the emission factors of various pollutants (Querol, 2016; Nielsen et al., 2017; Holder et al., 2019; Trubetskaya et al., 2021). The advantages of these studies are the controlled experimental conditions, the information about the combustion parameters (fuel type, combustion temperature, excess of air), and the opportunity to change these parameters, and thus EFs concerning a wide spectrum of fuels and combustion conditions were reported. However, since only a limited number of stoves and combustion scenarios were studied, it is difficult to extrapolate these results to a “real-world situation” of a city.

For this reason, other papers focus on the real-world situation and report the atmospheric concentrations of the BB-related air pollution. However, since the contribution of the traffic emission always interferes, the pure BB-related air pollution is difficult to study. Consequently, some studies selected specific locations like Glojek et al. (2022) in Loški Potok, Slovenia, that can be considered a model village of biomass burning emission with a negligible contribution of other sources of air pollution. Other studies utilize the source apportionment model (Aethalometer model) of Sandradewi et al. (2008) and reported BB- and FF-related BC concentrations separately (see, e.g., Dumka et al., 2018; Deng et al., 2020; Liakakou et al., 2020; Mbengue et al., 2020; Milinković et al., 2021).

Despite the reliable source apportionment of BC by the Aethalometer model, the determination of the source-specific EFs in a real-world situation is still problematic due to the lack of the CO₂ source apportionment. However, inverse modeling can offer the opportunity to trace back the air pollution to their sources. For example, Olivares et al. (2008) applied inverse modeling to retrieve the traffic- and BB-related emission factors of NO_X, PM₁₀, BC, and particle number.

In this paper we aimed to determine the biomass burning and traffic-specific BC emission factors in the urban atmosphere during the heating season. We used the carbon-balance method that required the simultaneous source apportionment of BC and CO₂ concentrations. The BC source apportionment was performed by the Aethalometer model (AM), while the source apportionment of CO₂ was implemented by multilinear regression analysis (MLR). After the source apportionment of both components, the specific emission ratios (ERs) for BB and FF were determined and converted to EF values following the carbon-balance method. The measurements were taken during a 3-month-long monitoring campaign in Ljubljana, Slovenia, during winter 2019–2020. The atmospheric concentration of black carbon was monitored with simultaneous CO₂ measurement at three locations of the city with different emission characteristics involving traffic- and heating-related emissions as well as an urban background site.

In the following we introduce our combined Aethalometer model–multilinear regression analysis (AM-MLR) method that we applied for the determination of the source-specific emission factors. We present the BB- and FF-related emission factors for three different locations of the city. In order to validate the AM-MLR method, an auxiliary measurement campaign was performed during summer, when only fossil fuel combustion was assumed to be present. The FF-related emission factors determined during the summer campaign were compared to the result of the AM-MLR method.

2 Methods

2.1 Measurement sites and instrumentation

The measurement campaign took place from December 2019 to March 2020. Three measurement sites were selected in the city with different microenvironments, source profiles, and emission activities (Fig. 1). One location was selected in the historical center of the city (Trnovo, TRO), where wood combustion represents the primary energy source for domestic heating during winter. This site is located in the restricted traffic area of the old town, where low direct vehicle emission is expected. The measurement setup was installed in a family house, with the sampling inlet on the roof, 8 m above the ground, to ensure the access of air masses arriving from all directions.

Another location was selected close to major roads and far from biomass burning sources that ensured measurement of the higher relative contribution of traffic emission. The instruments were installed in a waterproof cabinet in the open recreational area of the Atlantis sport complex in the BTC commercial center (BTC site). Due to the open environment of the location, a 3 m sampling height was chosen.

The third measurement location was the atmospheric observatory of Aerosol d.o.o (Skylab, SKY) that is considered to be an urban background location. This location is far from the major roads of the city and is not affected directly either by traffic emission or wood combustion. The sampling inlet was at 10 m above the ground, which ensured free access of air masses from each direction.

The equivalent black carbon (eBC, referred to as BC in the following) concentrations were monitored using multi-wavelength Aethalometers (AE33, Magee Scientific/Aerosol d.o.o. Slovenia, Drinovec et al., 2015) that measure the light attenuation of the particle sample collected on a Teflon-coated glass filter tape (M8060) at seven wavelengths (370–950 nm). The absorption coefficient of the particle sample (b_abs, Mm⁻¹) was obtained by dividing the attenuation coefficient by the multiple scattering parameter (C=1.39 for the M8060 filter tape; Weingartner et al., 2003; Yus-Díez et al., 2021). The BC mass concentrations were generated as the ratio of the absorption and the wavelength-dependent mass absorption cross-section parameter (MAC_λ, m² g⁻¹) provided by the manufacturer. Although the default MAC values cannot be used universally, their validity has been proven for Ljubljana (Ogrizek et al., 2022).

Aerosol size selection was provided at the inlet of the Aethalometer by a cyclone sampling head with a PM_2.5 cut-off diameter. The flow rate was set to 5 L min⁻¹ and the measurement time resolution to 1 min. The “dual-spot” technology enables the real-time loading effect correction, which is especially important when the spectral dependence of optical absorption is used for source apportionment (Drinovec et al., 2015).

The CO₂ concentrations were measured by flow-through CO₂ sensors (Carbocap GMP 343, Vaisala, Finland). The CO₂ sensors were directly connected to the exhaust of the AE33, and thus they analyzed the identical air stream to the Aethalometer. The accuracy of the sensor was 3 ppm + 1 % of the reading below the 1000 ppm concentration range, which was the case during the campaign even on the most polluted days. The response time of the sensor was comparable with the AE33, so the 1 min average signals of BC and CO₂ were well correlated when common sources were measured.

The three measurement systems were compared in the air quality laboratory of Aerosol d.o.o before the campaign. The variation between the AE33 units was below 1 % at 1 min averaging time for all the used wavelengths (880 nm for the BC concentration and 470 and 950 nm for the source apportionment). This precision was expected from the results of Cuesta-Mosquera et al. (2021), who compared 23 AE33 units and found variation between the measurement results of less than 1 %. The unit-to-unit variability of the CO₂ sensors was below 4 ‰ on a 1 min time basis.

2.2 Meteorological situation

The measurement campaign started on 6 December 2019, during a warming-up period that was continued by an unusually warm and dry January and February (Fig. 2). Table 1 summarizes the basic climatological anomalies compared to the reference long-term averages of the 1981–2010 period. The average monthly temperatures in Ljubljana during the 3-month-long campaign were warmer than the long-term averages of the 1981–2010 period (2.3, 1.7, and 4.8 ^∘C above the long-term average in December, January, and February, respectively). February 2020 was the second-warmest February in the history of measurements. Usually, January and February are the driest periods of the year, and in 2020, they were even drier than the average. The snow cover was negligible, limited to a few days during the measurement period.

Atmospheric dilution and dispersion significantly affect the pollution accumulation in the planetary boundary layer and play an essential role in the formation of the concentration level (Alfoldy and Steib, 2011). For the quantification of the dispersion of the pollution, the ventilation coefficient (VC, m² s⁻¹) as the product of the horizontal wind speed and the depth of the planetary boundary layer was applied. These parameters were provided by the Real-time Environmental Applications and Display System (READY) of the NOAA. The meteorological data used for the model calculations were obtained by the Global Data Assimilation System (GDAS), with a spatial resolution of 1 ^∘ and a temporal resolution of 3 h.

The VC value for Ljubljana followed a lognormal distribution from 166 to 63 500 m²s⁻¹ during the measurement period, with the median of the 2530 and 1460–4570 m² s⁻¹ interquartile range (see Fig. 3).

We identified the well-mixed, diluted cases based on high VC values, while low VC values referred to the opposite cases, when the atmospheric dispersion of the pollution was moderate, and thus emissions of local sources dominated the measured concentrations.

Table 1 Monthly meteorological anomalies relative to reference long-term averages of the 1981–2010 period.

2.3 ER and EF calculation

Air pollution emission from combustion sources is usually reported with respect to burnt fuel mass and is given in a fuel-consumption-specific EF in g(kg fuel)⁻¹. Although the combustion is never complete, more than 99 % of the fuel carbon content is oxidized to carbon dioxide (EEA, 2019) that can be used as a tracer for fuel consumption estimation. Dividing the pollution concentration increment by the CO₂ increment of the plume, the pollution-to-CO₂ ER can be determined.

The ratio of concentration increments in the plume can be calculated in an integrative or derivative way. If the time resolution of the measurement technique of the two components differs significantly, the two concentrations would not be correlated even if they have a common source. In this case, the integrative method is the preferable option for ER calculation. In this way the time integrals of the concentration peaks are calculated (peak area), and the ratio of the net peak areas (after background removal) provides the ER (see, e.g., Ježek et al., 2015). The disadvantage of this method is that the peak identification is arbitrary, and the background definition and removal burden the calculation with an additional uncertainty.

On the other hand, if the time resolution of the two measurement techniques is similar, the recorded pollutant concentrations originating from a common source are correlated in time. In this case a threshold value can be defined for the minimum required R² of the correlation. Above the threshold R² the two components are considered to originate from the same source (plume event), and the slope of regression provides the ER (derivative way). The offset of the regression line depends on the background concentrations that do not need to be taken into consideration during the calculation. This method also provides a well-defined plume event identification since the correlation between the two components is weak out of the plume, while a strong correlation indicates simultaneous concentration peaks of the two components.

In our case the $\mathrm{BC}/\mathrm{CO}$₂ ER was calculated by the derivative method, and later it was transformed to EF using the carbon content of the concerned fuel:

$$\begin{array}{}\text{(1)}& \mathrm{EF}\left[\mathrm{g}{\left(\mathrm{kg}\mathrm{fuel}\right)}^{-\mathrm{1}}\right]=\mathrm{ER}\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}/\mathrm{ppm}\right)\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{1.82}}}\cdot {\displaystyle \frac{\mathrm{44}}{\mathrm{12}}}\cdot \mathrm{CC},\end{array}$$

where CC is the carbon content of the fuel that is 0.86 for diesel oil and petrol (Huss et al., 2013), while it ranges between 0.42 (eucalypt) and 0.47 (olive wood) when considering dry wood (Gonçalves et al., 2012). Here we applied the value of 0.45, which corresponds to the most usually used pine and oak wood. The measured CO₂ concentration was converted from ppm to mg m⁻³ using a 1.82 mg m⁻³ ppm⁻¹ conversion factor considering the AMCA (Air Movement and Control Association International Inc.) atmospheric standard (T=21.11 ^∘C, P=1013.25 mbar) that was also applied by the Aethalometer for the BC concentration calculation. Molecular weights of CO₂ (44) and C (12) were used to calculate the carbon mass fraction in CO₂.

Table 2 Statistical metrics of the measurements at the three monitoring locations. Mean (with the source-specific percental share with respect to the total BC), standard deviation (SD), their ratio (coefficient of variation, CV), the three quartiles (1Q, “Median”, 3Q), minimum and maximum values as well as their difference (“Range”) were calculated from the FF- and BB-related BC and CO₂ concentrations. Statistical values for the temperature (T) and relative humidity (RH) were given as well for the BTC locations only.

2.4 Source apportionment and source-specific emission ratios

Measurement of the spectrally resolved absorption coefficient provides an insight into the composition of light-absorbing particles, allowing us to distinguish the highly (and widely) absorbing black carbon (soot) particles from brown carbon (light-absorbing organic aerosols; see Bond and Bergstorm, 2006; Drinovec et al., 2015). Fossil fuel combustion generates mostly pure soot particles that are strong light absorbers over the whole NIR-visible wavelength domain, while particles generated by biomass burning contain other light-absorbing compounds such as brown carbon that have characteristic absorbance bands in the near-UV domain (Sandradewi et al., 2008; Helin et al., 2018).

Sandradewi et al. (2008) developed the so-called “Aethalometer model”, where the absorptions at the 470 and 950 nm wavelengths were expressed as the sum of the absorptions of the FF- and BB-related BC components (BC^FF and BC^BB), while the ratios of the absorptions at different wavelengths follow a reciprocal power law of the wavelength ratio with a corresponding exponent (absorption Ångström exponent, AAE) of FF- or BB-related BC. In this study, an AAE of 1.15 was used for the FF-related emission whose value was determined during the summer auxiliary measurements when only FF sources were considered (see Sect. 2.5). For the BB-related component, an AAE of 2.1 was set according to the maximal AAE values that we measured at the TRO location (wood burning site) during nights when a low traffic contribution was assumed. The solution of the equation system results in the BB-related absorption at the 950 nm wavelength whose ratio to the total absorption provides the ratio of the BB-related BC concentration.

2.4.1 CO₂ source apportionment

In order to apply the carbon-balance method for the source-specific EF calculation, source apportionment of the carbon dioxide is needed as well, which was implemented using the BC source apportionment combined with MLR. The method assumes that either the FF- or BB-related CO₂ component is correlated with the corresponding BC component (BC^FF or BC^BB) in the plume. The total measured CO₂ can be expressed as follows:

$$\begin{array}{}\text{(2)}& {\mathrm{CO}}_{\mathrm{2}}\left(t\right)={\mathrm{CO}}_{\mathrm{2}}^{\mathrm{FF}}\left(t\right)+{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{BB}}\left(t\right)+{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{bg}},\end{array}$$

where CO${}_{\mathrm{2}}^{\mathrm{FF}}\left(t\right)$ and CO${}_{\mathrm{2}}^{\mathrm{BB}}\left(t\right)$ stand for the FF- and BB-related CO₂ components of the plume, respectively, while CO${}_{\mathrm{2}}^{\mathrm{bg}}$ represents the background concentration that changes much more slowly than the combustion-related components; thus, it can be considered constant during a plume event.

Equation (2) can be formulated using the FF- and BB-related BC concentrations and emission ratios (ER^FF, ER^BB),

$$\begin{array}{}\text{(3)}& {\mathrm{CO}}_{\mathrm{2}}\left(t\right)={\displaystyle \frac{{\mathrm{BC}}^{\mathrm{FF}}\left(t\right)}{{\mathrm{ER}}^{\mathrm{FF}}}}+{\displaystyle \frac{{\mathrm{BC}}^{\mathrm{BB}}\left(t\right)}{{\mathrm{ER}}^{\mathrm{BB}}}}+{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{bg}},\end{array}$$

or written in an equivalent form:

$$\begin{array}{}\text{(4)}& {\mathrm{CO}}_{\mathrm{2}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{{\mathrm{ER}}^{\mathrm{FF}}}}\left[{\mathrm{BC}}^{\mathrm{FF}}\left(t\right)+{\displaystyle \frac{{\mathrm{ER}}^{\mathrm{FF}}}{{\mathrm{ER}}^{\mathrm{BB}}}}\cdot {\mathrm{BC}}^{\mathrm{BB}}\left(t\right)\right]+{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{bg}}.\end{array}$$

Equation (4) expresses that the linear combination of BC^FF(t) and BC^BB(t) is correlated with the measured CO₂ using an appropriate ER^FF/ER^BB ratio. Our task is to find a particular ER^FF/ER^BB ratio which provides the best correlation between the two sides of Eq. (4). After the best correlation was found, the slope of the regression line provides 1/ER^FF, so ER^BB can also be calculated. The background CO₂ concentration determines the offset of the regression and does not need to be taken into consideration during the calculation. However, the background CO₂ provided by MLR is also valuable information that we present in this paper.

It has to be noted that oil burning for heating purposes is not usual in Ljubljana, so we could apportion all the FF-related BC to traffic sources. A high contribution of oil burning in the household energy production would interfere with the source apportionment that limits the applicability of the method in those locations, where oil heating is negligible.

The MLR analysis is a well-known and widely used method in source apportionment calculations; however, its combination with the Aethalometer model just recently appeared in the literature. Blanco-Alegre et al. (2022) applied the MLR method to decouple the BB- and coal-combustion-related BC based on source-specific correlations between specific tracers (K for wood burning and As for coal combustion).

Kalogridis et al. (2018) used the source apportionment information provided by the Aethalometer model for the source apportionment of carbon monoxide (CO) in Athens. They compared their result with the linear CO–NO_X model (see there) and concluded that the CO–NO_X model overestimates the BB-related CO contribution, maybe due to the photochemical loss of NO_X, while the MLR analysis provided more reliable results.

Table 3 Fitting parameters of the lognormal ER distributions (ng m⁻³ ppm⁻¹) at the three locations of the city (see Fig. 5). The distributions are normalized to 1. The derived median and mean ER values are also shown.

The AM-MLR presented here can thus be a universal technique for source apportionment of any air pollution component that co-emitted with BC (for example, organic carbon, CO₂, CO, NO, NO₂, SO₂, PM, or volatile organic compounds (VOCs).

For the application of the MLR analysis, the R statistical package “stats” (R Core Team, 2021) was used. The correlations were studied in a running time window with 1 h duration. During this time interval the background concentration is supposed to be constant, while the FF and BB sources have characteristic emission peaks.

The minimum R² criterion for the MLR analysis was set to 0.9, which represents a good correlation between the source-specific BC and CO₂ components. The correlation coefficient exceeded this threshold value during pronounced peak events only. A low R² value means uncorrelated BC and CO₂ peaks (i.e., shifted in time) or no presence of peaks in the time window. However, numerous BC peaks had to be discarded from the analysis due to the low or noisy CO₂ peak that resulted in a lower correlation coefficient than the threshold (see more in Sect. 3.2).

It should be noted that, during the 1 h time window, several FF- and BB-related sources contribute to the measured plume with different ERs. The emitted BC and CO₂ concentrations have been averaged out during the MLR, so the ER received from the actual time window refers to the 1 h average emission of the sources. The shorter the time window, the shorter the averaging period, which results in higher variation and a wider distribution of the ER values. However, the choice of the time window does not affect the mode of the distribution (the most frequent ER value).

In the following special conditions, the MLR method provided false results, and so they were discarded.

1. If the FF and BB components are well correlated (R²>0.8), the MLR method cannot separate the two components and provided similar ERs for the two components. Typically, this was the case when a transported pollution plume was measured within the FF and BB components arriving together at the measurement location, resulting in correlated concentration increments. In this case the ERs refer to the average BC emission ratio including all the combustion sources (FF + BB) and must be discarded from the results.

2. If the concentrations are dominated by one of the sources (BB or FF), good correlation was obtained between the corresponding BC component and the total CO₂ concentration. In this case the CO₂ source apportionment fails, and the total CO₂ increment is accounted for the dominant source, and consequently the calculation provides an underestimated EF. For this reason, cases when one of the components correlated well with the total CO₂ concentration (R²>0.8) were discarded from the analysis.

3. The maximal P value for significance criteria was set to 10⁻⁵ for both components. Results exceeding this threshold were discarded from the dataset.

Table 4 Fitting parameters of the lognormal EF distributions g(kg fuel)⁻¹ measured at the three locations under low (VC < 3200 m² s⁻¹) and high (VC > 4600 m² s⁻¹) ventilation conditions (see Fig. 6). The total number of cases (A) as well as the derived median and mean values were also presented. Bold numbers denote the EFs of direct sources, while italic numbers refer to EFs from transported pollution.

Table 5 Fitting parameters of the lognormal distributions of the fossil-fuel-related emission factor (EF, g(kg fuel)⁻¹) at BTC (all data) and the highway. The distributions are normalized to 1. The derived median and mean are also shown.

2.5 Auxiliary measurements

For the validation of the AM-MLR method, a well-defined case is needed with exclusively one type of source (traffic or wood burning). Since this was never the case during winter, we performed additional measurements during summertime next to the E61 highway ring around Ljubljana, where the plumes were expected to originate from pure FF emission sources only. A portable monitoring unit was used for the measurement, including an AE43 Aethalometer (Aerosol d.o.o, Slovenia) and a Vaisala GMP 343 CO₂ sensor, as in the winter campaign. The AE43 is a recently released battery-powered portable version of the AE33 Aethalometer with an identical optical chamber, flow system, and operation principle. In addition to its portable setup, the AE43 has a developed firmware and software system that offers improved user experiences with the real-time concentration and pollution-rose plots.

The measurement station was installed on an overpass road above the highway. The overpass makes a connection between two sections of an unpaved road that has negligible traffic (mostly agricultural vehicles), so practically the highway emission dominates the concentrations. Due to the fast fluctuation of the concentration and the short lifetime of the pollution peaks emitted by individual sources, a 1 s measurement time was used.

Since only FF-related sources were measured, the source apportionment and MLR procedures were not needed. The BC and the CO₂ concentration increments were well correlated during the peaks, and the slope of the regression was considered the ER^FF. Due to the rapid fluctuation of the concentrations, the regression was calculated using a 10 s running time window.

2.6 Uncertainty estimation

The resulting ER values are burdened by uncertainties that originated from several sources, such as (1) the concentration measurements, (2) the BC source apportion, and (3) the MLR calculation. In the following we discuss the different uncertainty sources and estimate the final error of the ER (EF) values.

The measurement is always burdened by a random fluctuation of the readings. As was already mentioned above, the unit-to-unit deviance was 1 % for the BC monitor and 0.4 % for the CO₂ sensor. In addition, there is a maximum of 2 % uncertainty of the CO₂ measurement and a significantly higher, 25 % systematic error of the BC measurement that comes from the uncertainty in the MAC value and the C factor (Ivančič et al., 2022).

Table 6 Emission factors of fossil fuel and biomass burning sources. Comparison of the results of this study with the literature. Results from the present study are shown as mean values and the interquartile range in brackets (q1–q3).

^∗ Converted from mg km⁻¹ units using the CO₂ EF from the same study. ${}^{\ast \ast }$ Converted from PM_2.5 g km⁻¹ EF using fuel consumption and BC percentage of PM_2.5 published by the same study.

The error of the BC source apportionment is caused by the uncertainties in the BB- and FF-related AAEs. Generally, when there is no independent measurement available for correlating the traffic or BB-related BC emissions, this error is estimated as 20 % (Healy et al., 2017). In our case, the source apportionment was combined with a correlation analysis (MLR) that increased the precision of the final values, since only those cases were considered where the source-specific BC–CO₂ correlation was high.

The abovementioned random error of the measured data and the uncertainty of the source apportionment result in uncertainties of the BC–CO₂ slopes given by the MLR analysis. Besides them, an additional error is generated by the failure of the presumed conditions of the MLR. During the MLR analysis it was assumed that the CO₂ background and the ER values are constant during the applied 1 h time window (see Eq. 4). Any significant variation of these parameters during the time window increases the uncertainty of the slopes.

The uncertainties of the slopes given by the MLR analysis performed by the “lm” function of the “stats” package were 6 % for the FF source and 8 % for the BB source on average. This is a combined uncertainty that includes the random error and the source apportionment uncertainty. It is seen that the original 20 % uncertainty was significantly improved by the MLR analysis due to the applied strict threshold of the correlation coefficient (R²=0.9). The systematic 25 % error of the BC measurement increases these uncertainties, so the final uncertainties of the ERs are 31 % for the FF sources and 33 % for the BB sources. The carbon content of the diesel fuel is exactly known, while the wood carbon content also has a ∼7 % uncertainty if the exact type of burnt wood is unknown (see Eq. 1 and the related text). So, the final EF^FF uncertainty is the same (31 %), while it is about 40 % for the EF^BB.

Table 7 Source apportionment of the combustion-related CO₂ and the background level (Bg) at the three monitoring locations as well as at the highway (FF component only). The mean (with the percental share with respect to the total combustion-related CO₂), standard deviation (SD), their ratio (coefficient of variation, “CV”), the three quartiles (1Q, “Median”, 3Q), minimum and maximum values, as well as their difference (“Range”) were calculated from hourly averages for the FF- and BB-related CO₂ concentration increments.

3 Results

3.1 Overview of the measurement results and diurnal cycle of the pollution

The statistical metrics of the hourly measurement averages at the three locations are summarized in Table 2. The BC^FF and BC^BB fractions are shown separately, as well as the CO₂ concentrations, temperature, and relative humidity. The meteorological parameters were measured at the BTC locations only but can be considered generally valid values for the whole city area.

It is seen that the traffic-related BC^FF component dominates the BC load at all the locations. The mean BC^FF concentrations were 3049, 2200, and 2650 ng m⁻³ at the BTC, SKY, and TRO locations, respectively, while the corresponding BC^BB concentrations were 1595, 1360, and 2180 ng m⁻³. The biggest difference between the FF- and BB-related components can be observed at the BTC location (64 % vs. 36 % of the total BC), while the smallest one was at TRO (55 % vs. 45 %), indicating a higher influence of wood combustion in the historical center of the city.

The spatial variation of the BC components shows an interesting pattern. Relative to the SKY location, the traffic-related FF component is higher by 38 % at BTC and 36 % at TRO. At the same time, the BB-related BC is higher by 17 % at BTC but 57 % at TRO, indicating that this (TRO) location is a definite hotspot in terms of wood combustion. On the other hand, the influence of traffic emission from the surrounding busy roads is still significant at the TRO measurement site even though it is located in a restricted traffic area.

The BC and CO₂ concentrations were highly affected by the atmospheric conditions at all the locations. Figure 3 shows the relationship between the VC and the BC concentration at the BTC and TRO locations. In the background of the figure the frequency distribution of the VC is plotted. The plotted BC values correspond to the average concentrations in the corresponding VC bin. The concentrations follow a decreasing trend with the increasing ventilation, indicating the dilution effect of the atmosphere. A significant concentration drop can be observed between 3200 and 4600 m² s⁻¹ VC values. It can be interpreted that, in the high concentration interval (BC > 4500 ng m⁻³, VC < 3200 m²s⁻¹), the local sources dominate the air pollution, while in the low concentration interval (BC < 2500 ng m⁻³, VC > 4600 m² s⁻¹) the contribution of transported, diluted pollution is the determinant.

The daily variation of pollution can be followed in the composite day concentration plots shown in Fig. 4. The FF- and BB-related BC concentrations are presented separately. It is seen that a pronounced FF peak can be found in the morning at 08:00 local time at all the locations, representing traffic emissions during the morning rush hours.

By contrast, the BB sources are more active in the afternoon. After 14:00 the BC^BB component starts to increase and reaches its daily maximum in the evening. An especially high evening maximum (3500 ng m⁻³) was found at the BB-influenced TRO location.

Figure 1 Measurement locations on the map of Ljubljana. HW shows the location of the traffic measurement next to the highway.

Figure 2 Time series of minimal, mean, and maximal daily temperatures, snow cover, daily sunshine duration and daily precipitation accumulation in Ljubljana from November 2019 to March 2020.

Figure 3 Frequency distribution of the ventilation coefficient (VC, histogram) and the averaged BC concentrations in the corresponding VC bins at the BTC and TRO locations (scatter plots).

Figure 4 Diurnal variation of the FF- and BB-related BC components as well as the CO₂ concentration at the (a) BTC, (b) SKY, and (c) TRO monitoring locations. The timescale represents local time.

3.2 $\mathrm{BC}/\mathrm{CO}$₂ emission ratios

Using the BC source apportionment results of the Aethalometer model, the MLR analysis provided the CO₂ source apportionment and the source-specific emission ratios. The normalized ER distributions are shown in Fig. 5 for the three locations. The distributions are wide and follow a lognormal pattern ranging from 10 to 1000 ng m⁻³ ppm⁻¹ according to the wide diversity of the sources. Lognormal curves were fitted on the distributions (solid lines in the figures), the parameters of which are summarized in Table 3. The mode and standard deviation that determine a normalized lognormal distribution are presented in the first two rows of the table. Since the median and mean differ from the mode for a lognormal distribution, these derived parameters are also shown in the last two rows of the table.

The wide distribution of ER can be explained by two main reasons. Firstly, the high variety of sources results in a wide range of emission ratios. For example, the BC emission factor of gasoline vehicles varies in the range of 0.001–0.01 g(kg fuel)⁻¹, while that of diesel vehicles falls into the 0.1–10 g(kg fuel)⁻¹ interval (EEA, 2019: 1.A.3.b.). Thus, the measured ER depends on the actual composition of the traffic, moving towards the higher values during the periods when the contribution of diesel sources (e.g., trucks, buses, and goods vehicles) is higher. In contrast, during periods when the traffic is dominated by personal vehicles, the ER decreases due to the higher contribution of gasoline vehicles.

Regarding the BB sources, the contribution of gas heating to the combustion-related CO₂ emission must be taken into account. The BC emission of gas heaters is much smaller than that of wood burning (0.6 vs. 74 g GJ⁻¹; EEA, 2019: 1.A.4.b), and thus the contribution of gas burning in the CO₂ plume dilutes the BB-related emissions. At the same time, the different burning conditions of wood stows from smoldering to high temperature flaming or the quality of the fuel (wood type, dryness degree) render high divergence of the emission ratios (see low fire–high fire variability in Table 6). More information about the relationship between combustion conditions and BC emissions can be found in the review of Shen et al. (2021).

Additionally, a measurement artifact caused by the high CO₂ background level also widens the ER distribution. Typically, the combustion-related CO₂ increments were measured in the 8–55 ppm interquartile interval, while the average CO₂ background concentration was 437 ppm with a 22 ppm interquartile range (see Table 7). This indicates how quickly a combustion-related CO₂ increment can immerse in the fluctuation of the background during the dispersion of the plume. For this reason, sources with high ER (i.e., low CO₂ increment) can be detected close to the sources only, and their relative contribution decreases with the increasing distance between the source and the measurement point. Therefore, diluted plumes always provide lower ERs than the direct ones, even if the composition of the sources is similar. Simultaneous measurement of direct and diluted plumes thus results in a wider ER distribution with a lower mode compared to the direct measurement. The same phenomenon leads to lower ER values in well-mixed atmospheres (high VC value) due to the dispersion of the CO₂ emission, while atmospheric inversion (low VC value) favors the detection of low CO₂ increments, thus resulting in higher ERs.

It is seen in the table that the ERs significantly vary between the locations. Maximal mean ERs were obtained at the BTC location (219 and 161 ng m⁻³ ppm⁻¹ for FF and BB, respectively), while the minimal mean ERs were found at the SKY location (160 and 100 ng m⁻³ ppm⁻¹ for FF and BB, respectively).

3.3 Emission factors of biomass burning and fossil fuel combustion

The emission factors were calculated from the ER values for biomass burning and traffic using Eq. (1). However, the dependence of the ER distribution on the pollution dispersion affects the calculated EF distributions as well. Figure 6 shows the EF distributions at the three locations considering two dispersion cases for both components, such as the (1) low-ventilation case (VC < 3200 m² s⁻¹) and the (2) high-ventilation case (VC > 4600 m² s⁻¹). The lognormal distribution function was fitted on the data points, and the fitting parameters are summarized in Table 4.

Figure 6a demonstrates that the atmospheric dispersion does not affect the FF emission factor distribution at the BTC site, while the mean BB emission factor is significantly lower in the case of high-ventilation conditions. This means that close FF sources were measured at the BTC site (local traffic), so the ventilation condition does not affect the EF. Regarding the BB component, a mixture of local and distant sources was measured, which later have more contribution with lower EFs during high-ventilation cases that shifted the EF distribution towards the lower EFs.

At the SKY location (Fig. 6b), no difference can be observed in the EF distributions, which indicates that diluted, distant plumes were measured even in low-ventilation cases (negligible contribution of local sources).

At the TRO location (Fig. 6c), both BB and FF emission factors decreased during high-ventilation conditions, but a more pronounced shift can be observed for the FF distribution.

The EF^BB mode under high-ventilation conditions at the BTC and TRO locations and the EF^FF mode at the TRO location are smaller than that obtained under low-ventilation conditions and match with the corresponding modes at the background (SKY) location (marked by italic numbers; see Table 4). This can be interpreted by the fact that transported plumes dominated the measurements during high-ventilation conditions at all locations and that a direct FF plume was measured at BTC only. It means that EF values under low-ventilation conditions (marked by bold numbers) are less affected by the atmospheric dilution and better estimations of the real emission factors.

Figure 5 FF-related (orange) and BB-related (blue) ERs at (a) BTC, (b) SKY, and (c) TRO locations, respectively. The distributions are normalized to 1. Lognormal distributions (solid lines) were fitted to the results: the parameters are summarized in Table 3.

Figure 6 FF-related (brown bars) and BB-related (blue bars) emission factors at the (a) BTC, (b) SKY, and (c) TRO locations, respectively. Low (VC < 3200 m² s⁻¹) and high (VC > 4600 m² s⁻¹) ventilation cases were plotted separately. Lognormal distributions (solid and dashed lines) were fitted to the results: the parameters are summarized in Table 4.

Figure 7 Distributions of the emission factors originating from traffic at BTC (all data) and the highway. Lognormal functions were fitted on the points with parameters summarized in Table 5. Colorized scatter symbols refer to the literature data (only X axis concerned). See more details in Table 6.

Figure 8 Distribution of the transformed FF EFs referring to the diesel emission at the BTC location (a). Scatter points show average EF data from the literature and the highway measurement (only X axis concerned). (b) Distribution of the BB EFs at the TRO location. Scatter points show average EF data from the literature and the BTC result (only X axis concerned). See more details in Table 6.

In Table 6 we summarize our results together with other data from the relevant literature. Mean values are shown in the table according to the literature data with the interquartile ranges in brackets. The presented data correspond to the low-ventilation case when the EF values were less affected by the atmospheric dilution. Notwithstanding, the comparison with the literature data is still problematic. Our results represent the average case of numerous urban sources involving low BC emitters (or non-smoking sources) that mostly contribute to the CO₂ increment (e.g., gas heating, gasoline vehicles). Thus, our results show lower EFs than of individual sources published in the literature. In the following we discuss our results in the context of the literature data considering the abovementioned aspect.

3.3.1 Traffic emission

Since the traffic-related EF does not depend on the ventilation condition at the BTC site, all the measured data were used for EF distribution without the consideration of the VC value. Figure 7 shows the EF distribution at BTC and at the highway that was measured during the summer campaign. The lognormal fits on the measured data are also shown. It is seen that the EF distribution at the highway site is much wider according to the applied short averaging window (10 s) during the correlation analysis that allows us to detect even individual sources. On the other hand, the two distributions covered each other with similar modes (0.33 and 0.36 g(kg fuel)⁻¹ at BTC and the highway, respectively; see Table 5). This good agreement between the AM-MLR method and the pure FF measurement verifies the validity of the AM-MLR method and indicates that the EF values were not distorted by the dilution effect.

In Fig. 7 relevant data from the literature are also shown in scatter plots (see more details in Table 6). Enroth et al. (2016) studied EFs of a mixed fleet in Finland near a highway. Their mean EFs were in the 0.15–0.54 g(kg fuel)⁻¹ range that overlaps with the EF distribution curve at BTC provided by the AM-MLR method.

Blanco-Alegre et al. (2020) measured BC EF in a 1 km-long urban tunnel in Braga, Portugal. Tunnels ensure well-defined conditions for traffic EF measurements with concentrated pollution that mostly originates from vehicle emission. The authors obtained an average EF of 0.31 g(kg fuel)⁻¹ for the fleet of nearly 56 000 vehicles, whose composition is probably similar to the Slovenian fleet (Cooper, 2020). This value is in very good agreement with the result of our AM-MLR method at BTC (0.39 g(kg fuel)⁻¹ average).

Olivares et al. (2008) measured source-specific black carbon concentration in Temuco, Chile, by Aethalometer and particle soot absorption photometer (PSAP). They determined the EF for a mixed fleet by inverse modeling that gave results of 0.35 g(kg fuel)⁻¹ mean EF for Aethalometer and 0.61–0.73 g(kg fuel)⁻¹ EFs for PSAP that fit into our EF distribution.

Assuming that BC emissions of gasoline vehicles are negligible compared to those of diesel engines (as is supported by tailpipe emission measurements – EEA, 2019), all the measured BC^FF can be attributed to diesel emission. On the other hand, the diesel-emission-related carbon dioxide can be estimated based on the share of diesel cars in the vehicle fleet, that is, 36 % in Slovenia (National interoperability framework – portal NIO, https://nio.gov.si/nio, last access: January 2023). This means that the diesel-emission-related CO₂ is roughly 36 % of the total CO${}_{\mathrm{2}}^{\mathrm{FF}}$. The emission factor of diesel engines can thus be calculated by dividing the original EF by 0.36.

In Table 6 the transformed EFs are presented for BTC, TRO and highway locations. These numbers refer to the diesel EF only, and they are in good agreement with Brimblecombe et al. (2015), who reported 1.28 g(kg fuel)⁻¹ diesel EF from a tunnel experiment in Hong Kong. The reported EF values from individual diesel cars (Ježek et al., 2015; Alves et al., 2015; Zavala et al., 2017; EEA, 2019) and individual truck emission monitoring (Ban-Weiss et al., 2009; Dallmann et al., 2011) are in good agreement with our transformed EF distribution (Fig. 8a).

3.3.2 Biomass burning

According to the literature data, the biomass burning EF from individual stove emission measurements ranges from 0.063 g(kg fuel)⁻¹ (Sun et al., 2018) to 0.83 g(kg fuel)⁻¹ (Holder et al., 2019; Akagi et al., 2011). The wide dispersion of the literature values indicates the high variety of BB EFs according to the stove type and combustion conditions. Figure 8b demonstrates that most of the literature data fall above our EF distribution measured at the TRO location. The lower EFs we found here can be the consequences of the contribution of gas combustion sources that are common all around the city. Gas burning emits a very small mass of aerosol particles compared to wood combustion, but at the same time, it significantly contributes to the CO₂ emissions from domestic heating. Since gas combustion for heating probably has the same time pattern as wood combustion (i.e., concentration increments during the evening and cold weather but drops during the midday and warmer periods), the CO₂ increments that correlate with the BC^BB component partially originated from gas heating. Our method thus cannot uniquely identify EFs from pure wood combustion but instead refers to the emission factor of the general domestic heating, including non-smoking sources as well. In an ideal case, when the measured sources were exclusively fueled by wood, the heating-related EF would equal the EF^BB; otherwise, the higher the contribution of gas heating, the lower the EF.

However, we also note that the real-world EF data published by Olivares et al. (2008) and the stove emission EF for pine wood by Sun et al. (2018) fall on the low end of the EF distribution measured at the TRO location.

3.4 Source apportionment of CO₂ emission

Using the source apportionment of BC and the calculated BC ER values, the BB and FF source-related CO₂ components can be retrieved. By subtracting the total combustion-related CO₂ increment from the measured CO₂ level, the non-combustion-related CO₂ level can also be determined.

Table 7 summarizes the statistical metrics of the BB and FF source-related CO₂ concentrations as well as the background level at the three measurement locations. In addition to the absolute mean values of the BB- and FF-related CO_2, their relative contributions to the total combustion-related CO₂ concentration are also shown as percentiles.

It is seen that the average background CO₂ concentration was the same (∼436 ppm) at all the locations. On the other hand, the source apportionment of the combustion-related CO₂ shows significant variation according to the environmental conditions of the locations. At the BTC location the FF-related CO₂ component is slightly lower than the BB component (47 vs. 53 %), while at the TRO location, the BB emission dominates the CO₂ level (62 %).

4 Conclusions

Atmospheric concentrations of black carbon and CO₂ were monitored in real time at three urban locations in Ljubljana, Slovenia, which had different impacts of traffic and wood burning during the winter heating season. The source-specific BC concentrations from the Aethalometer model were used to apportion the combustion-related CO₂ by coupling a multilinear regression method. The analysis presumed two combustion-related sources, namely, domestic heating (biomass burning) and traffic (fossil fuel combustion). The combined AM-MLR method provided consistent and realistic “real-world” emission ratios and emission factors for the three measurement locations. The method can be further generalized for source apportionment of other combustion-related components, and that of EFs can be determined later. Information about the source-specific EFs helps to estimate the pollution emission rates based on the fuel consumption.

The specific conclusions are as follows.

1. The traffic-related BC^FF concentration was higher than BC^BB at all locations. The smallest difference was found at TRO (wood combustion site), while the largest difference was obtained at BTC (traffic site). In contrast, the heating-related CO₂ concentrations were higher at all the locations.

2. The determined ERs follow a wide lognormal distribution according to the variety of the fuels (from the non-smoking gasoline or natural gas to the BC-producing diesel oil and wood) and sources (DPF equipped vs. conventional diesel vehicles; different types and conditions of wood stoves), as well as combustion conditions (high temperature, excess of air vs. low temperature, deficit air conditions). Also, it was shown that the distances of the sources affect the ER, since the relative contributions of high ER sources (means of low relative CO₂ emission) are lower for higher distances due to the dilution and fast dispersion of the related CO₂ increment.

3. Using the literature data of the carbon content of the fuels (diesel oil vs. wood), the related emission factors (EFs) were determined. The determined mean traffic-related EF (0.36 and 0.39 g(kg fuel)⁻¹ for TRO and BTC, respectively) is in good agreement with published EF values for a mixed traffic fleet. Using the relative ratio of gasoline and diesel fleet for Slovenia, the diesel-emission-related EF could be calculated (1.00 and 1.08 g(kg fuel)⁻¹ for TRO and BTC, respectively), which is in good agreement with diesel emission factors published in the literature.

4. The BB-related mean EF (0.13 and 0.16 g(kg fuel)⁻¹ for TRO and BTC, respectively) is lower than the majority of the relevant literature data reported for individual stoves. The difference might be explained by the CO₂ contribution of other, non-smoking combustion sources (i.e., gas heating).

5. The AM-MLR method was validated by direct traffic emission monitoring next to the highway during summertime, when only traffic-related sources were most likely sampled. Thus, the FF-related emission factors could be directly determined without source apportionment. The similarity of the modes of the two distributions (0.33 and 0.36 g(kg fuel)⁻¹ for BTC and HW, respectively) indicates that the AM-MLR method provided reliable results.