Firewood residential heating – local versus remote influence on the aerosol burden

We report the first-time use of the Lagrangian particle dispersion model (LPDM) FLEXPART to simulate isotope ratios of the biomass burning tracer levoglucosan. Here, we combine the model results with observed levoglucosan concentrations and δ¹³C to assess the contribution of local vs. remote emissions from firewood domestic heating to the particulate matter sampled during the cold season at two measurements stations of the Environmental Agency of North Rhine-Westphalia, Germany.

For the investigated samples, the simulations indicate that the largest part of the sampled aerosol is 1 to 2 d old and thus originates from local to regional sources. Consequently, ageing, also limited by the reduced photochemical activity in the dark cold season, has a minor influence on the observed levoglucosan concentration and δ¹³C. The retro plume ages agree well with those derived from observed δ¹³C (the “isotopic” ages), demonstrating that the limitation of backwards calculations to 7 d for this study does not introduce any significant bias. A linear regression analysis applied to the experimental levoglucosan δ¹³C vs. the inverse concentration confirms the young age of aerosol. The high variability in the observed δ¹³C implies that the local levoglucosan emissions are characterized by different isotopic ratios in the range of −26.3 ‰ to −21.3 ‰. These values are in good agreement with previous studies on levoglucosan source-specific isotopic composition in biomass burning aerosol. Comparison between measured and estimated levoglucosan concentrations suggests that emissions are underestimated by a factor of 2 on average. These findings demonstrate that the aerosol burden from home heating in residential areas is not of remote origin. In this work we show that combining Lagrangian modelling with isotope ratios is valuable to obtain additional insight into source apportionment. Error analysis shows that the largest source of uncertainty is limited information on isotope ratios of levoglucosan emissions. Based on the observed low extent of photochemical processing during the cold season, levoglucosan can be used under similar conditions as a conservative tracer without introducing substantial bias.

1 Introduction

Organic aerosol (OA) has anthropogenic and biogenic sources, being either released as primary OA (POA) or formed as secondary organic aerosol (SOA). Most of the anthropogenic emissions originate from combustion of fossil fuels or biomass. The biogenic particles are predominately SOA formed by the photo-oxidation of biogenic volatile organic compounds (VOCs). In the atmosphere, OA undergoes various physical and chemical processes, such as ageing by photolysis and photo-oxidation or deposition by sedimentation and wash-out. Particles have a direct radiative effect by absorbing and scattering solar radiation. Moreover, they act as cloud condensation nuclei (CCN), leading to cloud formation, which indirectly impacts the radiation budget. Being exposed to OA containing hazardous components, humans experience severe health impairments such as cardiovascular and respiratory diseases (Li et al., 2008, and references therein). Thus, OA affects air quality, health and climate.

Biomass burning is an important source of OA. Hallquist et al. (2009) estimated that biomass burning releases 42 Tg C a⁻¹ into the atmosphere, which is about a quarter of the global emitted particulate carbon. Such estimates are associated with considerable uncertainties. Parts of the uncertainties result from the lacking knowledge on source distribution and strength, as well as from the incomplete understanding of the loss processes. Since biomass burning substantially contributes to the OA hazards, it is of great scientific and societal interest to accurately apportion its sources and quantify its sinks. For the source apportionment of biomass burning aerosol, factor analyses, chemical mass balance and Lagrangian techniques are employed (e.g. Busby et al., 2016; Zheng et al., 2002). Chemical mass balance modelling used levoglucosan as the specific non-reactive molecular marker of biomass burning in aerosol (e.g. Fine et al., 2002) because it is only formed by the thermal breakdown of cellulose, and it is then emitted in large quantities. The accuracy of such studies is limited by considerable uncertainties in the emission factors and by the fact that levoglucosan was recently proven to be chemically unstable. Recent laboratory studies have shown that levoglucosan reacts with OH radicals within a lifetime of a few days under typical atmospheric conditions (Hennigan et al., 2010; Sang et al., 2016). This new finding opens up new potential applications, especially in the field of isotopic analyses.

Sources of biomass burning aerosol can be significantly better constrained by taking into account the stable carbon isotope ratio of levoglucosan (Gensch et al., 2014). This option is based on the fact that at the emission, levoglucosan has a source-specific isotopic composition, the “isotopic fingerprint”. Furthermore, chemical processing leads to isotopic fractionation due to the kinetic isotope effect (KIE), which is distinct for each reaction (making an “isotopic footprint”). Consequently, highly innovative source apportionment methods aim to combine trajectory and wind-based models with isotopic analyses, which deliver additional information for validation. Recently, Gensch et al. (2018) used isotopic measurements together with the Lagrangian particle dispersion model (LPDM) FLEXPART (Stohl et al., 2010) to investigate photochemical ageing processes in biomass burning aerosol. To this end, the photochemical age of particulate levoglucosan was derived from observed isotopic ratios, employing the “isotopic hydrocarbon clock” equation on the one hand and back trajectory simulations on the other. For the latter, a post-modelling numerical approach was developed to describe the mixing with freshly emitted levoglucosan and to quantify this impact on the isotopic composition. The results of these two independent methods agreed well on average. Moreover, the agreement between “isotopic” age and “retro plume” age demonstrates that the modelling results are not significantly influenced by limiting these to a few days for this study. Isotopes are thus a useful tool to evaluate the effect of the finite retro plumes. As a consequence, Gensch et al. (2018) showed that the degree of photo-oxidative ageing of particulate levoglucosan can be quantified by combining laboratory KIE studies and observed isotopic composition at sources and in the field, as well as back trajectory analyses. Yet, the scatter in the individual data pairs (model vs. observation) pointed out the need to improve the identification and distinction of contributions from different source types, which is possible by using the full dispersed output of FLEXPART.

As a particular form of biomass burning, home heating with firewood is a major contributor to the fine dust in the cold season in the mid-latitudes and high latitudes. According to the German Environmental Agency (UBA), small wood stoves in the residential sector provide only 1.5 % of the total energy supply but contribute 16 % to the total PM_2.5 emissions in Germany. This is comparable to the total road traffic PM_2.5 exhaust (Amann et al., 2018). For pollution mitigation, an accurate apportionment of local emissions versus remote transport is necessary. The main objective of this study was to implement stable carbon isotopes in the full dispersed output of FLEXPART by explicitly tracking the levoglucosan fraction containing ¹³C. In order to determine the model performance for given conditions, the sensitivity of the simulation responses to uncertainties of the governing atmospheric processes described in FLEXPART was examined. Finally, the set of selected modelling routines were applied in a case study with the goal to assess the contribution of local vs. remote emissions from firewood domestic heating to the particulate matter (PM) sampled at two measurement stations of the North Rhine-Westphalia Environmental Agency, LANUV. Thereby, the measured levoglucosan concentration and isotopic composition in the sampled aerosol were used to evaluate the model performance.

2 Experimental

The aerosol PM_2.5 fraction was sampled on quartz filters at two of the numerous LANUV monitoring network stations (Pfeffer et al., 2013). The sampling time was 24 h, and filters were changed daily at 00:00 UTC+1. For this study, two sampling sites with contrasting characteristics were chosen. The “urban background” station, hereinafter referred to as STYR, is situated in Mülheim-Styrum (51.453459^∘ N, 6.865050^∘ E; with site information available at https://www.lanuv.nrw.de/luqs/ messorte/pdf/STYR.pdf, last access: 27 September 2017), while the “remote/rural” one, hereinafter referred to as EIFE, is located in the hilly Eifel region (50.653234^∘ N, 6.281008^∘ E; with site information available at https://www.lanuv.nrw.de/luqs/messorte/pdf/EIFE.pdf, last access: 27 September 2017). For laboratory and model analyses, 25 pairs of aerosol filters collected on the same day in the cold seasons 2015–2017 at each of the two sites were decided on. The main criterion in selecting those was to provide a broad geographical coverage for the wind directions for the sampled air masses around the measurement sites. A list of the selected samples can be found in the Supplement (Table S1 in Sect. S1, containing the sampling dates, levoglucosan loading of the filters and the main origin of the sampled air masses estimated by HYSPLIT trajectory analyses).

Levoglucosan concentration was measured by ion chromatography (Kuepper et al., 2018) at the LANUV. Isotopic analyses were carried out at IEK-8, Forschungszentrum Jülich, by liquid extraction–thermal desorption–two-dimensional gas chromatography, coupled with isotope ratio mass spectroscopy (LE–TD–2D-GC-IRMS) (Gensch et al., 2018). Details on the experimental approaches are given in the Supplement (Sect. S2). Basic statistical analysis of the measurement results can be found in the Supplement (Table S3 in Sect. S3).

3 Modelling stable carbon isotopes with FLEXPART

Stable carbon isotope ratios of a VOC can be calculated in numerical atmospheric models by considering ¹²C and ¹³C isotopologues as separate species. These are treated individually during the simulations, taking into account the slightly different physical and chemical behaviour. Mostly, molecules containing ¹³C react a little slower than the light isotopologue due to the kinetic isotope effect (KIE). Denoting the rate constants for these two reactions by ¹³k and ¹²k, respectively, KIE is defined as ${}^{\mathrm{12}}k{/}^{\mathrm{13}}k$ and is typically expressed using the epsilon notation in parts per thousand (‰):

$$\begin{array}{}\text{(1)}& \mathit{\epsilon }=\left(\mathrm{1}-\text{KIE}\right)\times \mathrm{1000}\mathrm{‰}=\left(\mathrm{1}-{\displaystyle \frac{{}^{\mathrm{12}}k}{{}^{\mathrm{13}}k}}\right)\times \mathrm{1000}\mathrm{‰}.\end{array}$$

KIE can be experimentally determined in laboratory studies (Anderson et al., 2003; Sang et al., 2016). Carbon isotope effects are generally so small that the relation between the change in the isotopic ratios and the extent of chemical processing can be linearized without introducing significant bias (Rudolph and Czuba, 2000). Based on the proposed isotopic hydrocarbon clock equation, the photochemical age of a VOC can be determined from its isotopic composition when the source signature and the KIE of the atmospheric degradation reaction are known.

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_t={\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{0}}+t_{\text{av}}[\text{OH}{]}_{\text{av}}{k}_{\text{OH}}^{\text{OH}}\mathit{\epsilon },\end{array}$$

where δ¹³C₀ and δ¹³C_t represent the isotopic composition at the source and observation site, respectively, t_av[OH]_av is the average photochemical age, [OH]_av is the average OH concentration during the photochemical processing, k_OH is the rate coefficient of the species of interest with OH and ^OHε is the KIE of the latter oxidation reaction.

Not only are the fractionation effects small, but also the ratio of the rare isotope to the abundant one is very low. Therefore, carbon isotope ratios ${}^{\mathrm{13}}R=\frac{{}^{\mathrm{13}}\text{C}}{{}^{\mathrm{12}}\text{C}}$ are given using the delta value:

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}^{\mathrm{13}}\text{C}={\displaystyle \frac{{(}^{\mathrm{13}}\text{C}{/}^{\mathrm{12}}\text{C}{)}_{\text{sample}}-{(}^{\mathrm{13}}\text{C}{/}^{\mathrm{12}}\text{C}{)}_{\text{VPDB}}}{{(}^{\mathrm{13}}\text{C}{/}^{\mathrm{12}}\text{C}{)}_{\text{VPDB}}}}\times \mathrm{1000}\mathrm{‰},\end{array}$$

where ${(}^{\mathrm{13}}\text{C}{/}^{\mathrm{12}}\text{C}{)}_{\text{VPDB}}$ is the internationally accepted Vienna Pee Dee Belemnite (VPDB) value of 0.0111828 (Brand et al., 2010; Craig, 1957). The source-specific carbon isotope ratios of atmospheric trace organic components are introduced in the simulations considering the emission data for the investigated VOC (details in the following sections).

In this study, detailed information on origin and pathway of the collected aerosol particles was obtained by calculating retro plumes backwards from the sampling sites with the LPDM FLEXPART, version 10.2beta (Seibert and Frank, 2004; Stohl et al., 2010; source code available at https://git.nilu.no/flexpart/flexpart, last access: 4 October 2018). European Centre for Medium-Range Weather Forecasts (ECMWF) 3-hourly data with a resolution of ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$ on 91 vertical levels were used as driving meteorology (Owens and Hewson, 2018). For every investigated day, ca. 200 000 model particles were released hourly for 24 h at the measurement stations. Levoglucosan (LG) was implemented as an aerosol biomass burning tracer which is subjected to photochemical degradation by OH as well as to wet and dry deposition. The aerosol particle population defined in the input has a log-normal size distribution with $\stackrel{\mathrm{‾}}{d}=\mathrm{0.25}$ µm, ${\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{d}}=\mathrm{1.5}$ and a density of 1.4 kg m⁻³ (Fiebig et al., 2003). To simulate the wet aerosol particle removal, the new deposition module from Grythe et al. (2017) was activated, using three-dimensional cloud water fields from the ECMWF data. For below-cloud scavenging, a coalescence probability of 1 was set for both rain and snow (Grythe et al., 2017). For the chemical loss, the OH decay rate constant was set to $\mathrm{2.67}\times {\mathrm{10}}^{-\mathrm{12}}$ cm³ molec.⁻¹ s⁻¹ (Sang et al., 2016). The model was run backwards over 7 d. This is at the higher end of the expected levoglucosan lifetime for a mean OH concentration of 0.5×10⁶ molec. cm⁻³ in the cold season in Europe (Gerasopoulos et al., 2012; Rohrer and Berresheim, 2006). The output resolution was set to $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ on five vertical levels. The resulting retro plumes entail the source–receptor relationship reflecting the deposition- and decay-corrected receptor sensitivity to potential upwind sources. See Supplement (Sect. S4) for details.

To investigate the contribution of relevant domestic heating sources to the biomass burning aerosol sampled at the receptors, the retro plumes were folded with monthly mean gridded levoglucosan emissions during the cold season in Europe (Seibert and Frank, 2004). The emission inventories were derived from firewood consumption in the targeted European countries, population density and levoglucosan emission factors of firewood burned in common wood stoves (Akagi et al., 2011; Fine et al., 2004; Jimenez et al., 2017; Schauer et al., 2001) (firewood consumption by UN data: https://data.un.org, last access: 10 March 2017; population density data by NASA: https://neo.sci.gsfc.nasa.gov, last access: 6 February 2017). Details are given in the Supplement (Sect. S5). Domestic heating is the main source of the sampled levoglucosan in this study, since FIRMS (Fire Information for Resource Management System a NASA product: https://firms.modaps.eosdis.nasa.gov, last access: 10 January 2019) fire inventories show no larger scale open fires affecting the sampling during the considered periods. To determine the receptor sensitivity to home heating emissions, a dynamic footprint layer from 100 to 300 m height was considered (Hüser et al., 2017). The resulting folded retro plumes quantify the contribution of each individual source in kilograms per cubic metre (kg m⁻³) to the levoglucosan sampled at the measurement site. The corresponding concentration is then derived by summing up all contributions.

Levoglucosan δ¹³C at the sampling sites was determined by introducing ¹³C-LG as an additional model tracer (Gensch et al., 2011; Stein and Rudolph, 2007). It has the same physical and chemical properties as ¹²C-LG, except for a reduced OH reactivity due to the kinetic isotope effect. In the simulations, the rate of the ¹³C-LG chemical loss was derived by using KIE $=\mathrm{1.00229}\pm \mathrm{0.00018}$ (Sang et al., 2016). ¹³C-LG at the source was calculated using ¹²C-LG emissions together with the specified source isotopic ratios. The isotopic composition at the sampling point was derived as the ratio between the slightly different ¹³C-LG/¹²C-LG retro plumes folded with the corresponding ¹³C-LG/¹²C-LG emission inventories (for details, see Supplement, Sects. S4, S5).

4 Results and discussion

4.1 FLEXPART sensitivity studies

Different FLEXPART modules, describing chemical decay and dry and wet deposition, were successively activated. Modelling parameters, such as meteorology, levoglucosan lifetime and emission data, were varied to reveal the governing simulated processes and to assess the modelling performance. Following the overall goal of this study to separate local from remote sources of firewood domestic heating aerosol, changes in output depending on selected input parameter modifications were evaluated for consistency.

4.1.1 Driving meteorological input data

Unpredictability of the driving meteorology is one of the major error sources in Lagrangian modelling (Angevine et al., 2014; Davis and Dacre, 2009; Lin, 2013), leading to concentration uncertainties of up to 40 %. In this study, two global numerical weather prediction models, ECMWF and Global Forecast System (GFS, with the same horizontal resolution but with 39 vertical levels) (Global Climate and Weather Modeling Branch, 2003), delivered the input meteorological fields for two otherwise identical simulations with FLEXPART. The derived concentration and δ¹³C values were compared.

Tables S6.1 and S6.2 present model results obtained with the two meteorology sets (Supplement, Sect. S6). Differences in the isotopic composition are well below the experimental uncertainty of 0.6 ‰. Levoglucosan concentration calculations based on ECMWF agree fairly well with the measurements, showing mean deviations of 9.1 % and 4.2 % for the EIFE and STYR datasets, respectively. In general, the concentrations obtained when using the GFS meteorology are higher than for ECMWF initiated calculations, with a difference mean of 6.9 % and 4.6 %, respectively (Supplement, Figs. S6.1 and S6.2).

Figure 1 Simulated retro plumes (top; for the colour code, see legend) as well as centroid back-trajectory (bottom, solid lines) and the corresponding mixing-layer heights (bottom, dashed lines) for each release hour, using the ECMWF and GFS meteorology, exemplarily presented for the 18 February (left) and 5 March 2016 (right) samples. The FLEXPART analyses are done for the EIFE site (red star).

Two extreme cases are presented in Fig. 1. For 18 February 2016, the simulations show a 1.5 % deviation between the concentrations yielded using the ECMWF and GFS meteorology. This difference is largest (86.2 %) for 5 March 2016. In the former example, similar mixing heights (H_mix) and mixing events (when the centroid height H_traj<H_mix) are predicted. For the latter case, higher GFS H_mix are derived. Lower receptor sensitivity to emissions, expressed as t_res (Supplement, Sect. S5), is expected due to dilution. The t_res from ECMWF runs is smaller though. A possible explanation is that due to shallow mixing layers, a large part of the emitted particles likely leaves the investigated footprint layer and does not contribute to the receptor (Hueser et al., 2017).

Overall, the simulated concentration and δ¹³C based on ECMWF meteorological fields show good agreement with the values obtained when using GFS data, the uncertainties being in the same range as those of the experimental observations. Differences in the modelled levoglucosan concentration and isotopic composition are expected, mainly due to the different parameterization of sub-grid processes in the atmospheric boundary layer in the GFS and ERA5 models (Berrisford, 2011; Han, 2011). Mixing influences the concentration of a model tracer that is released in the planetary boundary layer. Yet, during the more stable cold season over Europe, the differences are expected to be small. Due to the higher vertical resolution and ability to more accurately account for topography, ECMWF meteorology was chosen to initialize the model for the further runs.

4.1.2 Photochemical ageing

Based on the FLEXPART full dispersed output of the “inert tracer” scenario, the retro plumes were divided into age classes, separating the data by the time the particles in the considered portion need until they reach the sampling point. Percentages of the total levoglucosan contribution to the sampled air were calculated. The results show that emissions during the last 24 h before sampling contribute on average 49 % of the sampled aerosol, being considered “1 d old”. Detailed model results are given in Table S7.1 (Supplement, Sect. S7). A total of 30 % of the contributed emissions occurred between 48 and 24 h and 10 % between 72 and 48 h before sampling, being considered “2 and 3 d old”, respectively. Only 11 % of the total emitted particles were older than 3 d. The simulations thus show that the major part of the sampled aerosol originates from local to regional sources being emitted during the sampling day and the day before. Exemplarily, Fig. 2 shows emission contributions of different ages for 1 April 2017. Due to the young age of the sampled aerosol and the typically low OH concentration in the cold season, photochemical ageing is not expected to be the governing loss process in this study.

Figure 2 Levoglucosan emission contribution to the EIFE site (green star), divided into age classes, for 1  April 2017. The plot shows source regions for 1 d old (35 % of the collected levoglucosan, red), 2 d old (35 %, yellow) and 3 d old (11 %, blue) particles, as well as particles older than that (20 %, grey).

To investigate the influence of the photochemical ageing of the sampling aerosol, simulations were carried out, implementing the chemical degradation of levoglucosan as described by Sang et al. (2016). Indeed, differences between the results of the inert and “reactive tracer” scenarios are on average 10 % and 7 % for the EIFE and STYR stations, respectively, both being within the experimental error range. Changes in the isotopic composition are at both stations 0.2 ‰ on average, being smaller than the measurement precision.

Overall, under the investigated conditions, both concentration and isotopic composition at the sampling site are rather determined by mixing with fresh emissions than by chemical loss processes. This agrees well with the study by Busby et al. (2016), which pointed out that levoglucosan is relatively stable during winter due to the low OH concentration.

4.1.3 Deposition

Simulations including dry as well as wet deposition were carried out. The calculated concentration and isotope ratios were compared with the reactive tracer (without deposition) scenario to quantify the contribution of dry and wet removal, respectively.

The initialized aerosol particle population has a mean diameter of 0.2 µm, with a corresponding average settling velocity of $\mathrm{5.6}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹. Since the aerosol in this study is relatively young, most of it stays in the accumulation mode. Thus, gravitational settling concentration losses are negligible, amounting to 1.5 % and 0.7 % at the EIFE and STYR stations, respectively. Consequently, changes in the isotopic composition are minor as well (details in the Supplement, Sect. S8).

Further, the existent simulations show that wet deposition removed minimal amounts of the emissions. This might be explained by a short exposure of aerosol to weak precipitation of less than 5 mm per 6 h in the investigated periods. Moreover, there were no in-cloud scavenging events for the investigated periods. Possible explanations are the absence of fog or spreading of emissions in layers lower than the cloud bottom layer height. Similarly to the low impact on concentration, wet deposition had no significant influence on the isotopic composition of the sampled aerosol either.

4.1.4 Domestic biomass burning emission estimates

Given the lack of levoglucosan emission data from residential heating, an approach was developed here to estimate these from available information. Therefore, population density, country-specific firewood consumption and levoglucosan emission factors for typical fuel used in residential heating were considered (details in the Supplement, Sect. S5). Accordingly, uncertainties arising from potential spatial and temporal variability of the emission intensities are brought into the calculations, such as errors due to seasonal and regional differences in the wood acquisition and consumption. By considering the monthly mean consumption of firewood, which is provided by a personal survey, weekly to diurnal variances in the emissions are neglected. Additionally, weighting with the population density does not reflect the real spatial distribution of the wood consumption. As an example, fireplace heating is rather unusual in cities with high population densities. The conversion factor from wood weight to levoglucosan emission depends on the wood type and the combustion process (Akagi et al., 2011; Fine et al., 2004; Jimenez et al., 2017; Schauer et al., 2001).

The injection height of a fire emission is usually parameterized based on exhaust magnitude and temperature. These are quite similar for domestic woodstoves, unlike in the case of large open fires. According to Zhang et al. (2014), the footprint layer for domestic heating emissions stretches from 100 to 300 m. The model results are not sensitive to the footprint layer height as long as it is inside an effectively mixed layer (Hüser et al., 2017). This is valid for the investigated period, since the FLEXPART-simulated mixing heights drop under 300 m in fewer than 10 % of cases. Furthermore, increasing the thickness of the footprint layer has no major influence on the model outcome because of two counteracting effects. A greater dilution reduces the impact of a source, while a wider spread of the emission increases the residence time of model particles in the footprint layer.

Figure 3 Distribution frequency of the observed levoglucosan concentration and δ¹³C.

Figure 4 Comparison between observed and simulated levoglucosan concentration at the EIFE and STYR sites. Also, the 1 : 1 line is given (dashed line).

4.2 Case study

Based on the sensitivity study findings, modelling routines and parameters were selected and applied in a case study aiming to assess the contribution of local vs. remote emissions from firewood domestic heating to the aerosol sampled at EIFE and STYR sites.

The measured levoglucosan concentration for the investigated samples varies over more than 1 order of magnitude (from ∼ 10 to 500 ng m⁻³), being lower overall at EIFE (54.2 ng m⁻³ on average) than at STYR (152.1 ng m⁻³ on average). The histograms depicted in Fig. 3 show further differences between the sites with a unimodal vs. a multimodal distribution of the concentrations observed at EIFE and STYR, respectively. According to this, there are two major categories of sources contributing to the sampled aerosol, i.e. regional and close-by sources. The urban site region is affected not only by regional upwind emissions, like the remote one, but also by sources very near to the receptor.

Measured versus simulated concentrations are depicted in Fig. 4. Generally, the model results are at the lower end of the observation range, showing mean absolute percentage deviations of 42.0 % and 53.8 % for the EIFE and STYR sites, respectively. For the rural station, the model overestimates the concentration of levoglucosan outside the experimental error ranges in 20 % of cases and underestimates it in 44 % of cases. By contrast, at the urban stations, there is only one case of model overestimation. Here, the model predominantly underestimates levoglucosan concentrations, far outside the experimental error ranges.

Figure 5 δ¹³C of the sampled levoglucosan at EIFE and STYR sites. The shaded areas represent ranges of observed levoglucosan source specific isotope ratios in aerosol formed during the combustion of soft (light green) and hard wood (dark green) (Sang et al., 2012).

Figure 5 shows the measured δ values of the sampled levoglucosan for both measurement stations. The isotopic composition ranges between −26.3 ‰ and −21.3 ‰. In theory, such high variance can be explained either by different source specific isotopic ratios of the contributing emissions or by a different extent of the chemical processing of levoglucosan in the sampled aerosol.

Simulations showed that most contributing emissions are only 1 d old. As a consequence, levoglucosan degradation is not expected to be significant for this study. This assumption is supported by the reduced photochemical activity in the dark winter months (Busby et al., 2016). Thus, mixing of sources characterized by manifold isotopic composition is the more likely explanation of the measured δ¹³C. The combination of isotope ratio with concentration ambient observations in form of a Keeling plot (see below) was employed to verify this conclusion derived from the model. The δ¹³C distribution shown in Fig. 3 is narrower for the EIFE site than for STYR, where mixing of individual different sources seems likely.

A model–observation comparison analysis gives the possibility to assign individual isotopic signatures to different source regions. Within this study, retro plume analyses providing the main wind direction were used to confine geographically source regions for the sampled aerosol (details in the Supplement, Sect. S9). Weighted mean measured δ values of these assigned source regions are given in Fig. 6. Within measurement error ranges, there is no significant difference in δ¹³C for the different source regions, indicating that there is on average no significant spatial dependence of the isotopic signature of levoglucosan emissions.

Figure 6 Isotopic signatures for different source regions, including corresponding sample number as well as the average δ¹³C. The source regions for the sampled aerosol were designated based on the main wind direction (see text).

4.3 Keeling plot analysis

To investigate the chemical stability of levoglucosan for this study, a Keeling plot approach was employed. The test hypothesis was that chemical decay plays an important role here. Accordingly, the Keeling plot describes the mixing of two reservoirs in this study of fresh, “isotopically light”, high-concentrated emissions with aged, “heavier”, low-concentrated background (Lin, 2013). For this analysis, the measured isotopic ratio was plotted vs. the inverse concentration (Fig. 7). A linear regression analysis was carried out. Remarkably, according to the 95 % confidence interval analysis, the yielded y-intercept range of −25.3 ‰ to −21.4 ‰, theoretically describing the sources, agrees well within error ranges with the published isotopic composition measured in aerosol from the combustion of various C₃ plants (Sang et al., 2012). A y intercept of $-\mathrm{23.2}\pm \mathrm{0.1}$ ‰ was derived, theoretically representing the mean isotope ratio of fresh emissions. The knowledge gained from the Keeling plot analysis was used to initialize the FLEXPART isotopic runs. Since the information on type of fuels consumed in certain regions is very scarce, a constant source-specific δ¹³C₀ of −23.2 ‰ was considered to calculate the ¹³C emissions used in the simulations (Sect. S5 in the Supplement). This isotopic ratio is in the range of soft woods (conifers). Since background levoglucosan concentration data were not available, the lowest measured concentration (12.4 ng m⁻³ at the EIFE station on 10 November 2015) was considered to be a constant background value. A corresponding δ¹³C of $-\mathrm{24.0}\pm \mathrm{0.3}$ ‰ was derived (see Fig. 7). The most important outcome was that the slope of the fitted line to the experimental data was found to be slightly negative, contradicting the assumption of an aged heavy background. This analysis thus demonstrated that the initial hypothesis of levoglucosan chemical instability was false. Consequently, the variability in the observed δ values is likely due to the contribution of local to regional sources that possess different isotopic ratios in the above-mentioned range and not due to chemical processing.

The isotopic ratio of the assumed background, lower than that of sources, cannot describe the “expected” photochemically aged aerosol. This result calls into question whether the initially postulated reservoirs explain the observed concentrations and delta values. Background levoglucosan can originate either from a diffuse source that is not related to the population density or from air masses that are not taken into account by the 7 d retro plumes. The former hypothesis is not plausible in winter. The latter is likely in terms of the concentrations but not compatible with the delta values, since the isotopic ratio of 7 d older levoglucosan increases by 2 ‰. The Keeling plot is therefore a strong indication that a significant background cannot be reconciled with the observed isotope values. Thus, a systematic underestimation of the source strength might be a better explanation for the observed data. In the first model runs, the background levoglucosan concentration and isotopic ratio were set to 12.4 ng m⁻³ and −24.0 ‰, respectively. Subsequently, model runs were repeated by gradually reducing the background. The best fit between the observed and calculated concentrations was reached by setting the background to zero (Supplement, Table S8.3). The linear regression analysis additionally shows a significant difference in the emission factors between EIFE and STYR. Nevertheless, emissions are underestimated by a factor of 2 on average.

Figure 7 Keeling plot depicting the observed levoglucosan δ¹³C vs. the inverse concentration (black symbols). Model results are given (red symbols), as well as the line fitted to the experimental data (dashed). The light and dark grey shaded areas represent the 95 % prediction and confidence intervals, respectively.

The modelled δ¹³C results shown in Fig. 7 are located next to the line fitted to the observations. Based on the isotopic hydrocarbon clock equation (Eq. 3), the isotopic age was calculated from modelled and observed levoglucosan isotopic ratios, using a δ¹³C₀ of −23.2 ‰ and kinetic data from Sang et al. (2016). The mean OH concentration was considered 0.5×10⁶ molec. cm⁻³. This yields a “negative” age for the observations (Supplement, Table S7.3), probably a consequence of the inaccuracy of the emission isotope ratio used, which is in the upper range of source-specific δ¹³C₀. Here it is to be noted that a levoglucosan lifetime that is 4 times shorter, when using the k_OH reported by Hennigan et al. (2010), would translate to even stronger emissions, characterized by isotopic ratios 2 ‰ to 3 ‰ lighter. Then again, these are at the lowermost end of the reported δ¹³C₀. The trajectory model has an age of around 1.5 d, with a small but very significant difference between EIFE and STYR (EIFE is older, as expected). The observed delta values show a similar difference under the premise used here that the emissions contributing to EIFE and STYR have the isotopic composition. This is to be expected for the trajectory analyses in this study, since the mean age is greater than 1 d and therefore is regionally averaged. Moreover, due to the number and duration of the sample collection, the averaging includes also a variety of prevailing conditions. In reverse, when assuming that the trajectory analyses correctly reflect the mean age, one can determine the emission isotopic ratios from the observations. Noticeably here, the values for EIFE and STYR are almost identical, supporting the good averaging conclusion.

Overall, this model–observation comparison study shows good agreement between the model and observations. This demonstrates that FLEXPART describes the atmospheric processes investigated in this study well. Despite unknowns expected to introduce biases in the analyses, FLEXPART simulations provide a good description of the sources and background.

5 Conclusions

In this study we have combined Lagrangian particle dispersion modelling with concentration and isotopic measurements of levoglucosan in PM sampled during the cold season at two LANUV stations as an innovative tool to investigate home heating aerosol sources and aerosol fate. To this end, we have successfully implemented ¹²levoglucosan and ¹³levoglucosan as separate chemical species in the LPDM FLEXPART, to calculate the isotopic ratio distribution of the specific biomass burning tracer from the source to the sampling site. The analysis of the full dispersed model output, in combination with emission inventories using the folded retro plume technique, yielded very detailed information on the source–receptor relationships. Thus, aerosol source contribution to the receptor sites and its loss processes during the atmospheric transport are quantified.

Sensitivity studies show that for this special case, varying model variables, such as meteorology, rate constant of the photo-oxidation reaction or deposition processes, does not yield significant changes in the simulation results. This is justified by similar vertical mixing parameterization in the wind models during the more stable cold season over Europe. Lower OH concentrations cause less photochemical degradation of levoglucosan. This can also be associated with the young age of the sampled aerosol derived from the simulations. Furthermore, sedimentation is an insignificant loss process for the fresh biomass burning aerosol. The activation of the wet deposition module leads to minor levoglucosan concentration reductions due to missing strong precipitation events in winter. Levoglucosan isotope ratios close to δ¹³C₀ of emissions can be either explained by deposition or by local to regional sources and dispersion. Since dry and wet deposition are insignificant, the latter hypothesis is the most likely.

The presented case study shows reasonable agreement between modelled and observed data within error ranges. However, the frequent underestimations, especially at STYR, might indicate unidentified sources or flaws in the levoglucosan emission strength. This comparison supports the fact that sources,which are very close to the receptor but not accurately described in the developed emission inventory approach strongly influence the local aerosol burden, particularly for the STYR site. The Keeling plot analysis proved that an aged background was insignificant for this study, supporting the conclusion that long-range transport minimally impacted the investigated aerosol. This can also explain the few overestimations of the derived concentration at EIFE, which might be caused by overestimated background levels rather than underestimated removal (Grythe et al., 2017). Repeated calculations reducing the background and increasing the emissions indicated that in the initial model runs, the source strength is most likely underestimated and the background levels overestimated. The measured δ¹³C values show much higher variability compared with the simulated isotopic ratios. This can be explained by possible individual source-to-source variation (e.g. due to differences in the fuel used). Within the variability of δ¹³C of emissions, the retro-plume-modelled age of levoglucosan agrees well with the age resulting from the observed isotopic ratios. This agreement demonstrates that the limitation to 7 d backward calculations does not create any significant bias. Finally, since observations as well as the retro plume analyse show that chemical ageing does not play a significant role in the cold season in central Europe, levoglucosan can be used as a “conservative” tracer under similar conditions. All these findings demonstrate the FLEXPART fitness to simulate aerosol processes occurring between the source and receptor. The sensitivity studies revealed individual factors leading to potential biases, while the comparison between simulated and observed concentration assessed the most probable sources and loss processes for the investigated aerosol.

Both sensitivity and case studies unquestionably point out that local/regional domestic heating is the major source contributing to the biomass burning aerosol burden under the investigated conditions. Thus, we have demonstrated that the developed modelling strategies are suitable to assess sources of biomass burning aerosol in living areas in winter. Under similar conditions, i.e. cold season in Europe, photochemical decay is negligible; therefore levoglucosan can be employed as an inert tracer in source–receptor studies, without introducing considerable bias.

For the future, more isotopic measurements of fuels used for domestic heating in Europe are essential to better constrain the isotopic signature of individual sources. Global modelling together with more frequent ambient measurements is necessary to describe the concentration and isotopic composition of background aerosol more accurately. Further studies, preferably on summer fires, are needed to test processes described by FLEXPART. Isotopic information will likely deliver the additional information to quantify aerosol photochemical ageing in a OH-radical-rich atmosphere, as well as the wet deposition during strong precipitation, leading to heavy levoglucosan removal.