Modelling anthropogenic aerosol sources and secondary organic aerosol formation: a wintertime study in central Europe

Anthropogenic aerosol particles remain a significant air quality concern in central Europe, particularly during winter months. This study employs the chemistry transport model COSMO–MUSCAT (COnsortium for Small scale MOdelling coupled with the MUltiScale Chemistry Aerosol Transport) to investigate particulate matter sources, with a focus on emissions from residential heating. The model results are compared with winter measurements from one site in Germany and two sites in the Czech Republic, where solid fuels are commonly used for heating. A non-reactive tagging method tracking primary organic matter (OM) reveals a high contribution from residential heating. Although the magnitude and temporal changes of the model results mostly agree with total OM values at two measuring stations, it appears to underestimate measurements at a site in the central Czech Republic. This underestimation is partly attributed to the inadequate representation of secondary organic aerosol (SOA) emitted from wood combustion. The study highlights the impact of anthropogenic volatile organic compounds (AVOCs) on SOA formation, which are currently underrepresented in air quality models. Sensitivity tests adjusting SOA yields and AVOC emissions increase OM concentrations by up to 40 % at the measurement sites. These findings emphasise the need for accurate parameterisation of AVOC-derived SOA formation and residential heating emissions to better tackle wintertime air quality challenges in central Europe.

1 Introduction

According to the European Environment Agency's air quality report, 238 000 premature deaths in the European Union (EU) in 2020 could be attributed to exposure to particulate matter with an aerodynamic diameter of 2.5 µm or smaller (PM_2.5) (EEA, 2022). In a review summarising multiple decades of research, Anderson et al. (2011) emphasise that exposure to PM significantly increases the risk of cardiovascular and respiratory diseases, posing a major global public health challenge. A report by the World Bank Group (2022) estimates that the societal cost of ambient fine particulate matter pollution in the European and central Asian regions reached 4.6 % of gross domestic product (GDP) in 2019. This estimate reflects the economic impact of PM_2.5-related health outcomes, including premature mortality, morbidity, and lost productivity. The target of the EU's Air Pollution Action Plan is a 55 % reduction in premature mortality due to PM_2.5 by 2030 compared to 2005 levels (EEA, 2022). However, based on self-reported data, 19 EU Member States still fell at least 30 % short of their 2030 PM_2.5 emission reduction targets in 2021. A study by Beloconi and Vounatsou (2023) found that as of 2021, 47.5 % of Europeans were living in areas where annual mean PM_2.5 concentrations exceeded the new EU limit of 10 µg m⁻³, which will come into force in 2030. Ground level measurements of PM_2.5 from the European Air Quality Monitoring Network for 2021 and 2022 show a striking gradient between clean and polluted areas. Eastern European regions and the Po Valley in Italy have the highest annual mean concentrations, while central and western Europe have much lower PM levels (EEA, 2019). Simulations of PM_2.5 exposure and PM_2.5-related mortality for the year 2015 by Gu et al. (2023) also indicate higher concentrations and associated health risks in eastern Europe.

In this transition zone between less and more polluted regions, the rural background station Melpitz in eastern Germany recorded the highest annual mean PM₁₀ concentration in 2021 as reported by the European Monitoring and Evaluation Programme (EMEP) (Fagerli et al., 2023). Previous studies in Germany have shown that long-range transport from eastern Europe, particularly from combustion processes, is a major contributor to regional background particle concentrations (van Pinxteren et al., 2019, 2016). The inflow of air masses from the east was associated with PM₁₀ concentration peaks, leading to an increase in exceedances of the current daily limit value of 50 µg m⁻³ (van Pinxteren et al., 2019). However, the relative contributions of multiple combustion sources to primary and secondary particles, as well as their transboundary transport, remain insufficiently quantified. This needs to be better characterised to enable effective and better targeted mitigation strategies to address the prevailing air quality challenges.

Source apportionment (SA) studies aim to link ambient concentrations of pollutants to their emission sources. Within chemical transport models (CTMs), two main methods can be used to do this: the emission reduction impact method and the mass transfer method. The emission reduction impact method, or brute-force approach, assesses how pollutant concentrations respond to specific emission changes (Thunis et al., 2019). An extreme case, the “zero-out” method, sets emissions from selected sources to zero and estimates their maximum possible impact on ambient concentrations. This approach helps to assess the potential impact of emission reductions on air quality (Clappier et al., 2017). Despite its conceptual simplicity, this method is computationally intensive, and the results are highly sensitive to the reference scenario chosen.

The mass transfer method, implemented in CTMs as the tagged species approach, estimates contributions from different source sectors and regions by tracing the mass transport of pollutants from emission sources to local concentrations (Thunis et al., 2019). In this method, new tracers are introduced for the pollutants of interest and labelled according to their emission sources, allowing them to be monitored throughout the model run (Kranenburg et al., 2013). This approach facilitates the study of source contributions across both spatial and temporal scales, with source definitions directly linked to the emission inventories used as model inputs (Mircea et al., 2020). Tagged species modules for particulate source apportionment are implemented in several chemistry transport models: e.g. LOTUS-EUROS (Kranenburg et al., 2013), DEHM (Brandt et al., 2013), Particulate Matter Source Apportionment Technology (PSAT) in CAMx (Yarwood et al., 2007), Tagged Species Source Apportionment (TSSA) (Wang et al., 2009), and the Integrated Source Apportionment Method (ISAM) (Kwok et al., 2013) in CMAQ (US EPA Office of Research and Development, 2024). Tagging approaches are not designed to assess the effectiveness of mitigation measures or the impact of emission reductions because they do not consider indirect chemical effects (Thunis et al., 2019). However, they do provide a direct, additive source attribution of pollutant mass concentrations.

The Transport and Transformation of Atmospheric Aerosol over Central Europe with an Emphasis on Anthropogenic Sources (TRACE) project aims to develop a comprehensive understanding of the contribution of transported anthropogenic aerosols relative to local emissions, integrating expertise in synergistic measurement methods and modelling tools. As part of this effort, this study aims to improve the understanding of the interaction between dispersion and transformation processes by investigating an area of large PM_2.5 concentration gradients in central Europe. Therefore, we implemented a non-reactive tagged species approach into the COSMO–MUSCAT model (COnsortium for Small scale MOdelling coupled with the MUltiScale Chemistry Aerosol Transport) (Wolke et al., 2012). The tagging approach is applied to identify primary PM sources with a focus on winter combustion emissions. Online and offline measurements from an extensive campaign in 2021 are used to validate the simulations and to improve the understanding of the local air quality.

2 Observations and modelling

2.1 Sampling sites

The TRACE winter campaign took place from 5 February to 24 March 2021 at three measurement sites in central Europe: two of the stations, Melpitz (DE) and Košetice (CZ), were already well established as part of the Aerosol, Clouds and Trace Gases Research Infrastructure (ACTRIS) and the European Monitoring and Evaluation Programme (EMEP), while the third (Frýdlant, CZ) was specifically installed for this project. The sites were selected to capture an important area of transition between polluted and less polluted regions in central Europe (see Fig. 1).

The research observatory Melpitz (51.54° N, 12.93° E; 86 m a.s.l.) is located 50 km east of Leipzig, Germany, observing atmospheric background conditions in central Europe. It has been operated by TROPOS for more than 30 years (Spindler et al., 2001; Poulain et al., 2011). The station is surrounded by grassland and flat agricultural land without any notable wind obstacles. About 60 % of the time throughout the year, the prevailing wind direction is south-west. These air masses are of maritime origin and reach Melpitz after having crossed western Europe and, in the immediate vicinity, the city of Leipzig. Easterly winds occur 17 % of the time throughout the year, bringing dry continental air masses affected by long-range transport from Poland, Belarus, Ukraine, Slovakia, and the Czech Republic (Spindler et al., 2001, 2012, 2013).

The National Atmospheric Observatory Košetice (49.35° N, 15.05° E; 534 m a.s.l.) is situated 60 km south-east of the Prague metropolitan area in the Czech Republic. There are several small settlements in the vicinity of the station; however, the district is one of the least populated in the country (Zíková and Ždímal, 2013). Surrounded mainly by agricultural land and some woodland, a medium-sized timber factory equipped with a biomass furnace is located 7.5 km from the site (Schwarz et al., 2016). In winter, air masses reaching the site predominantly originate from the south-west (44 %), passing over central Europe (Pokorná et al., 2022). Similarly, Lhotka et al. (2025) observed that the contribution of continental air masses was higher in winter compared to other seasons, highlighting a distinct seasonal difference.

The Frýdlant temporary measurement site (50.94° N, 15.07° E; 366 m a.s.l.) was set up 2 km north of the centre of Frýdlant, Czech Republic, close to the Polish border. The station is located on the north-western edge of the Jizera Mountains and is surrounded by forests and farmland. The Turów coal mine, a large Polish open-pit mine, is about 10 km south-west of the site. Lignite from the Turów mine is used to fuel the nearby Turów power station.

During the campaign period, the COVID-19 pandemic continued to affect Europe, with containment measures still in place. In Germany, there was a lockdown from 13 December 2020 to 3 March 2021. Non-essential businesses, schools, and childcare facilities were closed, and employees were required to work from home wherever possible. Essential services, such as supermarkets, pharmacies, and healthcare facilities, remained open. From 3 March 2021, restrictions were adjusted locally based on infection rates and other factors (BMG, 2023). In the Czech Republic, strict restrictions were in place until 11 April 2021 (Slabá, 2022). In Poland, a partial lockdown was enforced from 28 December 2020 to 14 February 2021. Some restrictions were eased on 1 January 2021, allowing shops in shopping centres and cultural institutions to reopen. However, on 20 March 2021, stricter measures were reintroduced until 9 April (A3M Global Monitoring GmbH, 2023).

2.2 Measurement data

A multi-device setup for data acquisition was in place at all three stations. The data presented in this study were measured with the instruments listed in Table 1. The instrumentation included an aerosol mass spectrometer for measuring the non-refractory near-PM₁ chemical composition (organic, nitrate, sulfate, ammonium, and non-sea-salt chloride) and a multi-wavelength Aethalometer for measuring the equivalent black carbon (eBC), connected to a dry PM₁₀ inlet. The mass concentration of PM_2.5 was measured by gravimetric filter sampling using a Digitel high-volume aerosol sampler with pre-heated quartz fibre filters. Samples were collected for 12 h, covering daytime (05:00 to 17:00 UTC) and nighttime (17:00 to 05:00 UTC). The filters were further analysed with a Sunset Laboratory thermal–optical transmittance (TOT) instrument according to the EUSAAR2 temperature protocol (Cavalli et al., 2010). Online Sunset Laboratory OC–EC data are also available for Frýdlant and Košetice. Carbon parallel-plate diffusion denuders were used to remove volatile organic compounds to prevent positive sampling artefacts caused by adsorption of gas-phase organics onto the filter (Turpin et al., 2000). An aerosol chemical speciation monitor (ACSM) was used for aerosol mass spectrometry in Košetice, while aerosol mass spectrometer (AMS) instruments were used in Melpitz and Frýdlant. Hereafter, we use AMS/ACSM to refer collectively to measurements from all three instruments deployed at the sites. The AMS/ACSM instruments measure total organic matter (OM), which we can compare directly with our model output, while the two Sunset Laboratory instruments detect organic carbon (OC). For better comparability, OC was converted to OM using an OM$/$OC ratio from the literature. Poulain et al. (2011) estimated an OM$/$OC ratio of 1.64 at the Melpitz station in winter 2009, with almost no diurnal variation. For a winter campaign in 2020 in Košetice, Pokorná et al. (2022) found a ratio of 1.51 (±0.36). In this study, we have applied a factor of 1.6 to the conversion of all Sunset Laboratory data. Polycyclic aromatic hydrocarbons (PAHs) were detected on the filters by gas chromatography–mass spectrometry (GC–MS).

Table 1 Measuring devices used at the three sites to obtain the data for this study.

2.3 Model description

All simulations were carried out with the multiscale model system COSMO–MUSCAT. It consists of two online coupled components, the regional numerical weather forecasting model COnsortium for Small scale MOdelling (COSMO) (version 5.05, Schättler et al., 2018) in conjunction with the aerosol chemistry transport model MUltiScale Chemistry Aerosol Transport (MUSCAT) (Wolke et al., 2012), developed at TROPOS. The model system is designed for aerosol chemistry process studies and air quality applications at the regional scale (Hinneburg et al., 2008; Heinold et al., 2011; Tõnisson et al., 2021; Wolke et al., 2012) and participated in model intercomparisons such as the Air Quality Model Evaluation International Initiative (AQMEII; Im et al., 2015; Galmarini et al., 2021). COSMO is a nonhydrostatic atmospheric numerical weather forecasting model based on the primitive thermo-hydrodynamic equations describing compressible flow in a moist atmosphere. The atmospheric equations are solved based on a terrain-following grid with rotated coordinates (Schättler et al., 2018). The meteorological model provides all the necessary meteorological fields (e.g. wind, relative humidity, temperature) to MUSCAT, which then simulates the transport and chemical transformations in the atmosphere for different gas- and particle-phase species. Transport processes include advection and turbulent diffusion, while physical loss processes are characterised by dry and wet deposition (Wolke et al., 2012). COSMO and MUSCAT operate largely independently on separate grids and are coupled at each horizontal advection time step (every 15–80 s), allowing highly time-resolved meteorological input for the chemistry transport model.

Anthropogenic emissions of atmospheric compounds are treated as prescribed point and gridded area sources. Emissions within Germany are provided by the GRETA database of the German Federal Environment Agency (UBA) (Schneider et al., 2016) for the year 2019 (resolution: 1 km × 1 km). For European emissions outside Germany, the CAMS-REG v5 emission inventory for the year 2018 (resolution: 6 km × 6 km) is used, provided by the Copernicus Atmosphere Monitoring Service (CAMS) (Kuenen et al., 2022). Emissions are treated according to the Gridded Nomenclature For Reporting (GNFR) (NFR-I, 2023); i.e. they are grouped into different emission sectors representing different source types (e.g. public power, traffic; see Table 2). The temporal variation in emissions (daily, weekly, and seasonal cycles) is accounted for by time profiles, which differ according to the emitting sector. These temporal profiles are largely based on those provided with the TNO_MACC-II inventory (Kuenen et al., 2014), with adjustments for livestock and agriculture emissions according to Skjøth et al. (2011). For GNFR sector C (“other combustion”), the temporal profile weight applied to the emission factor ranges from 0.37 to 2.54 over the study period (see Fig. A1 in the Appendix). Emissions are provided as aggregated totals for some pollutant groups, which we then break down into individual components. Primary particulate matter (PM) is split into elemental carbon (EC), primary organic matter (OM), sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), sodium, and other minerals. A further distinction is made between fine (< 2.5 µm) and coarse (2.5–10 µm) aerosol particles. Non-methane volatile organic compound (NMVOC) emissions are divided into 23 different hydrocarbon groups. The splitting profiles for PM and NMVOC are based on different literature sources and are also provided by CAMS (Kuenen et al., 2022). In this study, country-specific splitting profiles (based on the year 2018) are applied to the overall emission input (see Tables A2 and A1 for GNFR C splitting factors).

The emission of biogenic VOCs (BVOCs) is based on Steinbrecher et al. (2009) and improved for extended land use categories according to Luttkus et al. (2022). The primary natural aerosol components are emitted online in COSMO–MUSCAT. The estimation of desert dust mobilisation depends on soil texture and soil size distribution according to Tegen et al. (2002) and preferential source regions (Heinold et al., 2011; Schepanski et al., 2017) using the current wind fields and hydrological conditions provided by COSMO. The emission of sea spray aerosol is based on Barthel et al. (2019).

Natural fire emissions (e.g. EC, OM, and primary PM_2.5) are provided as point sources for the year 2021 by the Global Fire Assimilation System (GFAS) (Kaiser et al., 2012). These emissions are resolved into 24 h mean values with a specific injection height for each point source.

Dry deposition is modelled using the resistance approach described by Seinfeld and Pandis (2006). Aerosol particles and trace gases are also removed from the atmosphere through wet deposition, subdivided into in-cloud and below-cloud scavenging. Both processes are parameterised by size-dependent particle capture efficiencies and corresponding gas uptake coefficients (Simpson et al., 2012).

To describe the gas-phase chemistry, an extended version of the Regional Atmospheric Chemistry Mechanism, RACM-MIM2-ext (Karl et al., 2006; Stockwell et al., 1997; Karl et al., 2009), is used. The mass-based aerosol population is described using a hybrid bulk–bin scheme. It comprises 25 prognostic aerosol particle tracers, including primary PM_2.5 and PM₁₀, primary OM, EC, sulfate, nitrate, ammonium, and secondary organic aerosol (SOA), as well as six bins for sea salt and primary marine organic particles (diameter range: 0.01–10 µm) and five desert dust bins (0.2–48 µm).

Secondary inorganic aerosol is formed through reactions between ammonia and sulfuric or nitric acid, which are generated from the gaseous precursor species sulfur dioxide (SO₂) and nitrogen oxides (NO_x) (Hinneburg et al., 2008). The partitioning between the particle and gas phases depends on the ambient atmospheric temperature and humidity. The implementation of this particle–gas partitioning follows the equilibrium approach described by Galperin and Sofiev (1998), utilising the methods proposed by Mozurkewich (1993).

The formation of SOA is described by the module SORGAM (Schell et al., 2001), extended to include additional biogenic volatile organic compound (BVOC) precursors from isoprene, monoterpene, and sesquiterpene oxidation and highly oxygenated molecule (HOM) formation from all considered BVOCs (Luttkus et al., 2022). The module uses the two-product approach described by Odum et al. (1996), which splits each SOA product class – comprising reaction products from aromatic precursors, alkanes, alkenes, α-pinene, and limonene – into two pseudo-products. For each, the formation of low-volatility products and their gas–particle partitioning are simulated.

Low-volatility condensable products are formed through oxidation of organic precursor gases by the OH radical, the nitrate radical NO₃, and ozone. The amount produced is determined by a product species (i)-dependent stoichiometric coefficient (α_i) in the specific reaction of the chemical mechanism (Schell et al., 2001). The SOA mass resulting from gas–particle partitioning is then calculated using a partitioning coefficient K_om,i for each low-volatility product species following Pankow (1994). The partitioning coefficient depends on temperature and is influenced by the molecular weight and saturation vapour pressure of species i. Each pseudo-product consists of a gas-phase and particle-phase product with different α_i and K_om,i. All information regarding the precursor VOCs, SOA class names in both the gas and the particle phases, and the reactions and stoichiometric coefficients can be found in Schell et al. (2001) and in the supplement of Luttkus et al. (2022). The total SOA yield (Y) resulting from the two previous steps can be calculated according to Eq. (1), where M_o is the total available absorbing organic matter (Odum et al., 1996).

$$\begin{array}{}\text{(1)}& Y=\sum _{i=\mathrm{1}}^n{Y}_i=M_{\mathrm{o}}\sum _{i=\mathrm{1}}^n\left({\displaystyle \frac{{\mathit{\alpha }}_i{K}_{\mathrm{om},i}}{\mathrm{1}+K_{\mathrm{om},i}{M}_{\mathrm{o}}}}\right)\end{array}$$

Over a range of organic mass concentrations M_o, a precursor gas will have a range of aerosol yields Y. The relationship between yield and organic mass concentration can be determined through chamber measurements. To model this relationship, a curve is fitted by selecting the optimal values of α₁, α₂, K_om,1, and K_om,2 within the two-product framework. The sum of all particle-phase products considered gives the total SOA concentration.

2.4 Model setup

The domains for the COSMO–MUSCAT simulations were chosen to cover the three measurement sites. To reduce computational costs for the targeted horizontal resolution in the measurement region, the model is nested twice. The innermost domain TraceD1 covers 317 × 204 grid cells with a horizontal resolution of ∼2 km (see Fig. 1). The vertical resolution for COSMO in TraceD1 is 50 layers with a maximum height of 22 km, while MUSCAT uses only the lowermost 27 layers, i.e. up to ∼6 km. A common grid nesting approach is used for the inner domains. The results of the larger domains are used as lateral boundary conditions on the inner domains. The meteorological initial and boundary conditions for the European domain (N0) are provided by reanalysis data of the CAMS global atmospheric composition forecasts (Inness et al., 2019). The simulation covers the period from 1 January to 31 March 2021, including a 1-month spin-up, with an output resolution of 1 h. The model system is re-initialised every 48 h using the aerosol and trace gas concentrations at the end of the previous run and a 24 h COSMO pre-run to spin up the meteorology in order to avoid long-term drifts in the modelled meteorology.

Figure 1 Domains for COSMO–MUSCAT runs and localisation of the three rural background sampling sites (© OpenStreetMap contributors 2017, distributed under the Open Data Commons Open Database License (ODbL) v1.0) (OpenStreetMap contributors, 2017).

2.5 Source attribution in COSMO–MUSCAT

A source attribution module has been developed for COSMO–MUSCAT 5.05 to analyse the influence of specific source regions, point sources, and emission sectors on primary particulate matter compounds. This new tagging method allows for the individual tracking of emitted source-specific species during a single model simulation, thus enhancing the analytical capabilities of the model. Unlike the “zero-out” method, which requires multiple simulations for each source sector or region of interest, this new approach eliminates this need. As a result, the analysis is faster and less computationally intensive. To do so, an additional tracer is introduced into the model emissions for each species of interest from each defined source sector or source region and combinations of both (see Table A3 in the Appendix). This additional tracer is labelled with the source information and then processed in parallel. In this way, the concentration of each of the so-called tagged tracers is available in each grid cell of the model and at each time step. This provides detailed spatial and temporal information about the source contribution to local tracer concentrations. In addition to the concentration of each tagged species, the total concentration – representing the cumulative impact of all sources – is also available. This allows the relative contribution of each tagged source to be effectively calculated. An overview of the selected source sectors is given in Table 2. Tags for source regions can be specified via a text-based input file in which each surface grid cell can be assigned a region name. For this study, we have tagged emissions from all countries within the inner domain TraceD1, namely Germany, Poland, and the Czech Republic. Additionally, a source region “boundary” is introduced, referring to the transport from the coarser domains to the inner domain. Input from outside the European domain is not included in the “boundary” tagged sector.

Table 2 GNFR source categories considered in this study.

Transport (advection, diffusion, sedimentation) and removal (dry and wet deposition) processes are automatically applied to tagged tracers in the same way as for all other tracers. However, gas-phase chemistry and aerosol chemistry are not considered at present within the tagging algorithm. Therefore, only chemically passive tracers can be tagged, i.e. using the non-reactive tagging approach. This enables a high-spatial- and high-temporal-resolution analysis of the source composition of primary particles. As this study focuses on winter combustion processes, anthropogenic EC and OM emissions are tagged. EC and OM emissions are split into fine and coarse aerosol; therefore the same split is applied for the tagged tracers.

3 Results

3.1 Meteorology

During the campaign, notable meteorological events and sharp temperature changes occurred in a short period of time. In early February, a low-pressure system with cold air in the north and warm air in the south moved southwards, transporting cold air to the Balkans and Greece. On 7 and 8 February, strong easterly winds and heavy snowfall led to significant snow drifts in some areas of central Europe. This was followed by a week of clear nights with prevailing westerly winds and temperatures dropping to −20 °C. The model successfully captured the period of low temperatures and the snow event at all three stations (see Fig. 2). The snow event was followed by a cold episode, resulting in more stagnant conditions with a change in wind direction and decreased wind speed at all stations.

Figure 2 Meteorological parameters for the three stations. The top row shows surface temperature and precipitation, with the shaded area representing surface dust concentration in the size class <80 µm. The bottom row displays the modelled surface wind speed and direction.

In mid-February, a nearly stationary high-pressure system transported warm air from the Sahara into central Europe, driving a rapid temperature increase of up to 20 °C within a week. An omega blocking pattern over eastern Europe facilitated the inflow of dust that accompanied the warm air, allowing particles to travel as far north as Scandinavia (Hoshyaripour, 2021; Haarig et al., 2022). This event significantly affected all three stations, resulting in elevated surface dust concentrations of up to 50 µg m⁻³ (see Fig. 2). The unusually high, spring-like temperatures persisted until the end of February. Another significant dust event occurred on 3 March, originating from the Sahara and affecting central Europe. Although our model successfully simulated dust uptake, surface concentrations during this event remained lower than those observed in mid-February. In mid-March, a shift to westerly winds brought low-pressure systems accompanied by widespread precipitation over Germany. This was followed by an intrusion of cold air from the polar regions, resulting in sleet, snow, and gusty winds (DWD, 2021a, b). By the end of March, atmospheric blocking patterns established stable and dry conditions.

3.2 Measurements and model capability

To validate the model's performance, we compared the modelled PM_2.5 concentrations and their components with observational data. In our model, PM_2.5 comprises mineral dust, sea salt, organic matter (including primary organic carbon and secondary organic aerosol), elemental carbon, ammonium sulfate, and ammonium nitrate. Figure 3 presents the time series of the total PM_2.5 concentration. The model does not accurately represent the magnitude of the concentration peaks, especially in Košetice, where the average modelled values are almost 10 µg m⁻³ below the observed values. The normalised mean bias (NMB) reflects the systematic bias and indicates a strong underestimation of PM_2.5 by more than 40 % in Melpitz and Frýdlant and −57 % in Košetice. The Pearson correlation coefficient (R) provides an indication of the overall fit of the trend, ranging from −1 to 1. An absolute value of exactly 1 means that a linear equation perfectly describes the relationship between the model and the measurements. Among the sites, the PM_2.5 trend is most accurately captured in Košetice (R=0.61) and least accurately in Melpitz (R=0.25). The root mean squared error (RMSE) quantifies the error between measured and modelled surface-level mass concentrations. Overall, the model RMSE is high, with values of 14.26 µg m⁻³ for Melpitz, 13.85 µg m⁻³ for Košetice, and 10.92 µg m⁻³ for Frýdlant. Together with the NMB, the high RMSE indicates that the model tends to underestimate concentrations during periods of high concentration peaks, as the RMSE is particularly sensitive to outliers. All statistical parameters are presented in Table 3. Im et al. (2015) analysed the performance of multiple models in simulating PM_2.5 concentrations as part of the AQMEII model intercomparison project. They found that most models systematically underestimated PM_2.5 at rural stations, with biases ranging from −2 % to −60 %. The COSMO–MUSCAT model performed relatively well, showing a bias of −24.82 %. However, all models struggled to capture wintertime levels, underestimating concentrations by more than 50 % across all regions. During the first 2 weeks of February, the TRACE campaign revealed the largest discrepancies between observed and simulated PM_2.5 concentrations, with most other tracers also underestimated. Strong easterly winds until 8 February facilitated long-range pollutant transport to Melpitz and Frýdlant. The snow event on 7–8 February led to a decrease in PM_2.5 concentrations in Melpitz by approximately 10 µg m⁻³. In Frýdlant, a slight decrease of around 4 µg m⁻³ was observed after the event, while in Košetice, concentrations even increased by about 4 µg m⁻³, indicating limited overall washout effects. Concentrations rose again after the snow event, peaking around 10 February. The snow event was followed by a cold episode with stagnant conditions, reduced wind speeds, and a shift in wind direction, leading to pollutant accumulation. The model may underestimate residential emissions due to missing temperature dependencies and unaccounted-for COVID-19 lockdown effects. Increased heating activity due to unusually cold temperatures and limited mobility combined with stagnant meteorology could lead to the observed underestimation of PM_2.5. Considering only data from 15 February onwards improves model performance, reducing the RMSE to 9.64 µg m⁻³ for Melpitz, 12.44 µg m⁻³ for Košetice, and 7.63 µg m⁻³ for Frýdlant. Additionally, the overall trend is captured more effectively, with R increasing to 0.79 in Košetice, 0.65 in Frýdlant, and 0.61 in Melpitz.

Figure 3 Time series for PM_2.5 mass concentration for the three stations. Filter data compared with modelled primary and secondary aerosol mass concentration. The timestamp for the filter data corresponds to the time of filter collection.

Table 3 Time-averaged measured and modelled mass concentrations and the associated root mean squared error (RMSE), correlation coefficient (R), normalised mean bias (NMB), and fraction within a factor of 2 of the observations (FAC2) for the whole campaign period. The modelled data were adjusted to match the measurement intervals before statistical analysis. Online refers to in situ measurements and offline to filter sampling. A factor of 1.6 was applied to the OC measured by the Sunset Laboratory offline instrument.

To gain a better understanding of the remaining discrepancies between modelled and measured PM_2.5, we can evaluate the accuracy for each individual PM_2.5 component (see Fig. 4).

Figure 4 Box and whisker plots for observations and modelled data at corresponding times during the campaign period. Rectangular boxes display the first and third quartiles. The line within each box represents the median value. Outliers are excluded. For OM and EC, all data were averaged to 12 h. For comparability with AMS/ACSM measurements, a factor of 1.6 was applied to the OC measured by the Sunset Laboratory instruments.

3.2.1 Mineral dust

The Saharan dust outbreaks likely influenced the total PM_2.5 concentrations during the TRACE campaign. In the model, the February event brought high dust loads for several days and led to dust deposition at all three stations (see Fig. 2). Lidar measurements in Leipzig recorded pure dust conditions, but below 3 km height, aerosol from continental Europe was likely mixed into the Saharan dust plumes (Haarig et al., 2022). This event had a rather short travel time (less than 2 d) before reaching Leipzig. For the March event, the model also shows dust reaching the three stations, though the loads were not as high as during the second event. Observations by Haarig et al. (2022) detected mixed pollution–dust conditions after air masses were transported over Spain and France, reaching Leipzig after 3–4 d. It is possible that the model underestimated surface PM_2.5 during these events, potentially due to limitations in the model domain or insufficient vertical mixing to the surface.

3.2.2 Nitrate and sulfate

In Melpitz, the model performs well for sulfate, with a correlation coefficient of 0.71 and a small bias (NMB $=+\mathrm{10}$ %), while nitrate is overestimated (NMB $=+\mathrm{51}$ %), though its temporal variability is captured reasonably well (R=0.62) (see Fig. 4a and b). In Frýdlant, the model shows moderate correlations (R≈0.40–0.46) and biases (NMB $=+\mathrm{29}$ % for sulfate and +40 % for nitrate) and low agreement within a factor of 2 (FAC2 < 50 %). Košetice exhibits the weakest agreement, with low correlations (R=0.16 for nitrate, R=0.36 for sulfate) and underestimations of both species (NMB =−26 % for nitrate and −51 % for sulfate). These results are broadly in line with model performance criteria reported in the literature, e.g. NMB within ±45 % for sulfate and ±60 % for nitrate (Huang et al., 2021) or NMB within ±30 % and R>0.40 (Emery et al., 2017). This indicates that the model captures the general magnitude and temporal variability of secondary inorganic aerosol concentrations across the domain reasonably well, despite some site-specific discrepancies (Table 3). The AMS/ACSM may underestimate total sulfate and nitrate concentrations in winter, when particle growth shifts part of the mass beyond the PM₁ range (Poulain et al., 2020), though these species are generally predominantly found in PM₁ (Zhang et al., 2023). Given their relatively small contribution to total PM_2.5 at our sites, it is unlikely that secondary inorganic aerosols are responsible for the discrepancy between the predicted and measured PM_2.5 aerosol mass concentrations.

3.2.3 Elemental carbon

EC concentrations show overall good agreement with observations (see Fig. 4d). Our model aligns more closely with the Aethalometer data in Košetice (RMSE: 0.66 µg m⁻³, NMB: +30 %) and Frýdlant (RMSE: 0.93 µg m⁻³, NMB: −7 %) than in Melpitz, where it agrees well with the offline Sunset Laboratory measurements (RMSE: 0.37 µg m⁻³, NMB: −23 %). The discrepancy between Aethalometer and Sunset Laboratory measurements arises from the different carbon fractions they detect: the Aethalometer measures optically absorbing carbon (black carbon) in PM₁₀, while the Sunset Laboratory instruments quantify elemental carbon (see Fig. A3 in the Appendix). Although differences in particle size cut-offs must be considered when comparing observations and model results, Poulain et al. (2011) found that around 90 % of the mass of elemental black carbon (eBC) in PM₁₀ is contained within the PM₁ fraction. Comparing across these different size classes should therefore be reasonable. In Košetice and Frýdlant, our model slightly overestimates EC concentrations, with NMB values of +136 % and +144 %, respectively. For winter 2019, Aethalometer measurements reported 0.98 ± 0.76 µg m⁻³ BC in Košetice (Lhotka et al., 2025), while Pokorná et al. (2022) found 0.92 ± 0.77 µg m⁻³ for winter 2020. In comparison, our averaged model result for 2021 was 0.76 µg m⁻³. In Melpitz, literature data show significant variability in BC concentrations. Atabakhsh et al. (2023) reported a value of 1.38 µg m⁻³ converted to PM₁ using a multi-angle absorption photometer (MAAP) during winter 2016/2017. Later, van Pinxteren et al. (2023) observed a marked decrease to 0.5 ± 0.41 µg m⁻³ in winter 2018/2019, likely reflecting reduced emissions and meteorological influences.

3.2.4 Organic matter

The modelled OM that we refer to hereafter is the sum of the fine primary organic aerosol (OM in PM_2.5), the total SOA, and OM from outside the European simulation domain. Primary OM accounts for approximately half of the total OM, with mean contributions of 44 % in Melpitz and 48 % in Frýdlant and a slightly higher share of 57 % in Košetice (see Fig. A2 in the Appendix). Panel c in Fig. 4 compares all available OM values for our campaign period. Across all three stations, the comparison to the Sunset Laboratory data shows a systematic underestimation by the model, with large negative NMB values: −73 % in Melpitz, −79 % in Košetice, and −67 % in Frýdlant. RMSE values are also high for Melpitz and Frýdlant (4.95 and 5.18 µg m⁻³) but improve notably when considering only data from 15 February onwards, decreasing to 2.87 and 3.85 µg m⁻³, respectively.

The underestimation of these values by our model seems to have a large contribution to the total PM_2.5 underestimation. The discrepancy between Sunset Laboratory and AMS/ACSM observations may partly result from the different particle size ranges each instrument targets: Sunset Laboratory samples PM_2.5, while AMS/ACSM captures only PM₁. Nevertheless, since organic aerosol is predominantly found in the submicrometer size range throughout the year (Poulain et al., 2020), the comparison between modelled PM_2.5 OM and AMS/ACSM PM₁ OM remains meaningful. However, other factors contributing to the discrepancy cannot be ruled out. AMS/ACSM instruments are particularly well suited for capturing temporal variability due to their high time resolution and finer size cut-off (PM₁). The Sunset Laboratory instruments provide an estimate of the total carbonaceous mass and are useful for assessing the magnitude of concentrations. It uses the same filters as the gravimetric reference method, allowing a more direct comparison to total PM_2.5 mass and offering a more complete picture of the aerosol burden.

In Melpitz and Frýdlant, the model aligns reasonably well with AMS/ACSM observations, with RMSE values of 1.17 and 2.01 µg m⁻³ and NMBs of −8 % and +18 %, respectively. Correlation is also relatively strong in Melpitz (R=0.60) but lower in Frýdlant (R=0.19), where the model fails to capture diurnal variability. The model underestimates the OM concentrations by AMS/ACSM in Košetice (RMSE: 6.48 µg m⁻³; NMB: −74 %) and also does not fully reproduce the diurnal variations (R=0.39) (see Fig. A2 in the Appendix).

The correlation coefficient is good for Melpitz (0.60), indicating a good simulation of the overall trend as well. In contrast, for the Sunset Laboratory data, the correlation coefficient is only 0.24. The range of concentration of the Sunset Laboratory data is underestimated by the model (NMB $=-\mathrm{73}$ %), which indicates that the model might be missing OM at this site. For Frýdlant, the diurnal patterns are not well captured by the model (R in comparison with AMS = 0.19), while the overall range of OM concentrations is also underestimated (Sunset Laboratory offline NMB = −67 %). In Košetice, the AMS/ACSM and both Sunset Laboratory instruments give consistent results, while the modelled data are noticeably lower. The AMS/ACSM detects lower values than our model, and modelled diurnal patterns also do not match the observation (NMB = −74 %, R=0.39). In comparison to the Sunset Laboratory filter measurements, the model shows a similar underestimation (NMB = −79 %). For Košetice, the same concentration levels for PM₁ and PM_2.5 size class OM indicate a dominance of fine aerosols, while there is few coarse-mode organic aerosols. The correlation coefficient of the model concentrations against the AMS/ACSM measurement is lower (0.39) than that of the filter samples (0.63). If only the SOA components of the modelled OM are taken into account, the correlation coefficient compared to the AMS/ACSM for Košetice decreases further to 0.18. For Frýdlant and Melpitz, calculating the correlation coefficient using only SOA gives similar results to using the total modelled OM concentration. The reduced correlation in Košetice when isolating SOA implies that the model underestimates secondary aerosol at this site, thereby negatively affecting the overall correlation.

Previous Sunset Laboratory filter measurements taken in Melpitz in winter 2018/2019 found an average value of 3.2 ± 3.2 µg m⁻³ (van Pinxteren et al., 2023). AMS/ACSM data for winter 2019 also give comparable values of 2.08±1.6 µg m⁻³ (Poulain et al., 2011), while measurements with an ACSM in winter 2016/2017 show higher values of 6.21 µg m⁻³ (Atabakhsh et al., 2023). For Košetice, a good characterisation of the site is also given by various previous studies. AMS/ACSM measurements provide average values of 3.13 µg m⁻³ in winter 2019 (Lhotka et al., 2025) and 4.55±4.40 µg m⁻³ in winter 2020 (Pokorná et al., 2022). Mbengue et al. (2018) found an average OC concentration in PM_2.5 of 2.85±1.91 µg m⁻³ for the period 2013–2016. For our study period, we found Sunset Laboratory average filter values ranging from 5.06 µg m⁻³ in Melpitz to 7.74 µg m⁻³ in Košetice, exceeding typical values reported for previous years. This suggests a strong influence of meteorological conditions on the overall concentration levels.

The discrepancy between modelled and measured PM_2.5 concentrations does not appear to be primarily driven by deviations in elemental carbon, sulfate, or nitrate concentrations. The overall good agreement between modelled and observed EC values, with correlation coefficients up to 0.61 and low bias, indicates a reliable simulation of primary combustion aerosol emissions. The contribution of secondary inorganic aerosol to total PM_2.5 is limited, and the discussed modelling uncertainties are likely not the main reason for the underestimation of total PM_2.5. OM is significantly underestimated, especially in Košetice (NMB =−74 %, R=0.39), which explains a large part of the PM_2.5 model bias. We hypothesise that the underestimation of secondary organic aerosol is a major source of error in total PM_2.5 simulations. A spatial variation in the model's performance is apparent, with similar trends observed in Melpitz and Frýdlant, whereas Košetice exhibits distinct behaviour. The dominance of fine particles in OM, suggested by nearly identical concentrations in the PM₁ and PM_2.5 size fractions, points to elevated levels of secondary particles. Given the general agreement with previous campaigns, we can also conclude that our campaign, although characterised by unusual meteorological conditions during the first weeks, represents typical aerosol conditions at these sites. Therefore, the underestimation could indicate a general underrepresentation of SOA during winter in this area in COSMO–MUSCAT.

3.3 Source attribution for elemental carbon and primary organic matter

Since the model accurately reproduces EC concentrations, which represent a primary anthropogenic aerosol component, we conclude that anthropogenic sources are well represented in the model, enabling reliable identification of source contributions. Additionally, with approximately half of the total OM comprised of POA, we infer that overall source profiles can be effectively identified by analysing primary OM and EC using the non-reactive tagging approach. The results, shown as relative contributions to primary OM and EC for the cold and warm periods (Fig. 5), underline the importance of long-range transport of particles. The source sector “boundary” represents transported particles from the outer model domains into the innermost domain where tagging is applied. During the warm period, long-range transport accounts for about 38 % of both EC and OM in Melpitz, illustrating the significant influence of particles originating outside the domain. In Košetice, the contribution is 23.8 % for OM and 22.8 % for EC, while Frýdlant has the lowest influence, with 16.6 % for OM and 14.6 % for EC. The prevailing wind regime and the basin-like topography of the Czech Republic reduce the influence of long-range transport in Košetice compared to the other two stations. The “boundary” contribution to fine OM and EC is only slightly higher than for Frýdlant, which is located in the middle of the domain (see Fig. 8). Backward trajectory analyses (HYSPLIT; Stein et al., 2015) indicate that during the high-PM-peak event in early February, stationary meteorological conditions resulted in minimal air mass transport to all sites. This effect is particularly pronounced in Košetice, where strong local stagnation can be observed.

Figure 5 Relative source contributions to primary organic matter in PM_2.5 and elemental carbon in PM_2.5. Top: cold period (5–16 February 2021), bottom: warm period (16 February–23 March 2021)

Figure 6 Spatial distribution of the average absolute contribution of emissions from the source sector Public Power to the concentrations of primary organic matter in PM_2.5. Areas with many coal-fired power plants are highlighted.

The study region, spanning parts of Germany, Poland, and the Czech Republic, is characterised by a high density of active lignite mines (see Fig. 6). Lignite, a particularly emissions-intensive fuel, is the energy source for many large power plants in this area. Germany and Poland host the largest number of coal-fired power plants with the highest total capacities in the EU (Alves Dias et al., 2018). Emissions from power plants used for electricity and heat production are categorised under the source sector “public power”. Despite its proximity to areas with a high number of coal-fired power plants, the “public power” sector contributes only a small share to the overall concentration of primary OM and EC. Tagging results for this sector, split by country of origin, are shown in Fig. 7, indicating that the peaks in Frýdlant are predominantly driven by Polish emissions. The proximity of the Turów lignite power plant largely explains the observed peaks, especially during periods of low wind speeds. During other periods, emissions from German and Czech sources dominate. The influence of coal burning on air quality in Frýdlant is further amplified by its use in domestic heating. In 2017, 47.7 % of Polish households with individual heating relied on coal (Macuk, 2019).

Figure 7 Elemental carbon in PM_2.5 concentration, broken down by country of origin and source sector. Left: “public power”, right: “other combustion”. The wind barbs represent wind direction and speed, with 12 h averages modelled for each station at surface level.

Figure 8 Spatial distribution of the average absolute contribution of emissions from source sectors “other combustion” (top) and “boundary” (bottom) to the concentrations of primary organic matter in PM_2.5.

The sector “other combustion” includes combustion processes of private households, in particular domestic heating processes with all fuel types. This sector has the biggest contribution, with up to 76.3 % for EC and 72.6 % for primary OM in Frýdlant. Contributions to fine OM from the “other combustion” sector are the highest in the Czech Republic and in urban agglomerations in Poland and around Berlin, Germany (see Fig. 8, right panel). The main contributors to the concentrations observed at the stations are emissions originating within the country where the station is located. However, Melpitz stands out with the highest proportion of contributions from cross-border emissions. Atabakhsh et al. (2023) carried out a positive matrix factorisation (PMF) analysis over a period of 1 year in Melpitz. They found the highest coal combustion contribution to POA under the influence of easterly continental air masses. Furthermore, they found a temperature and RH dependence for the factor consisting of aged SOA and highly oxidised OA in winter, with the highest concentrations observed at temperatures below 0 °C and RH above 80 %. They concluded that increased precursor emissions due to higher heating activities and amplified aqueous-phase chemistry lead to increased SOA formation. This could suggest a potential additional underestimation of the SOA formation rate in early February, as strong easterly winds were observed, followed by a subsequent cold period.

4 Discussion

Chen et al. (2022) conducted a multi-year PMF source apportionment study across various locations in Europe. They identified a coal combustion factor of primary OA at only two sites: Melpitz (data collected in 2016/2017) and the urban location Kraków (data collected in 2018). The strong seasonal variations in this factor suggest that it originates from residential heating emissions. The study also examined Košetice, where no coal combustion factor was detected; however, biomass burning accounted for 15.5 % of the total OA in winter 2019. Lhotka et al. (2025) conducted a PMF study with data also collected in 2019 in Košetice. They identified a coal combustion factor with the highest contribution of 5 % to total OA in spring, while biomass combustion contributed the most in winter (12 % of total OA). Both factors showed similar diurnal cycles related to domestic heating, and a strong correlation between levoglucosan and the biomass combustion factor was observed in winter, indicating a high proportion of wood combustion. During a particularly cold period in January 2019, an increased contribution of coal was observed, probably due to its increased use in private households for heating, given its higher calorific value compared to wood. The results are also consistent with those of Horník et al. (2024), who performed a PMF study with samples collected during the TRACE campaign for water-soluble organic compounds using NMR. They found a high residential heating contribution with coal markers, indicating additional coal combustion in early February in Košetice.

Pokorná et al. (2018) analysed changes in PM_2.5 composition and sources from the 1990s to 2009/2010 in Košetice. During this period, the dominant sources shifted from lignite combustion by power plants and oil combustion to residential heating, mainly with coal and/or biomass. In the Czech Republic, only 5 % of total coal consumption in 2019 was used in the residential sector, as part of the “other combustion” source sector (IEA, 2021). Hovorka et al. (2015) conducted a receptor modelling study in a residential area 64 km north-east of Prague in winter 2013 and estimated that wood burning contributed 49 % to the mass of fine aerosol. They found high correlations between contributions from wood combustion and levoglucosan and suggested that wood combustion in local boilers is common in suburban areas in the Czech Republic.

The landscape surrounding Košetice is mainly agricultural with scattered woodland, and the only direct sources of pollution are local roads and domestic heating (Zikova and Zdimal, 2016). It is plausible to assume high rates of wood burning, given the proximity of the timber factory. Levoglucosan, an aerosol tracer that is associated with biomass burning, measured during the TRACE campaign, shows the highest mean concentrations in Košetice (0.32 µg m⁻³) and the lowest in Melpitz (0.15 µg m⁻³). Polycyclic aromatic hydrocarbons (PAHs) are also good tracers of combustion processes; e.g. retene is a unique marker of wood combustion (Ramdahl, 1983). The average retene concentration in Košetice is 2.13 ng m⁻³, with an average total PAH concentration of 24.43 ng m⁻³. In Frýdlant, the average total PAH concentration is comparable (24.73 ng m⁻³), but retene concentrations are lower (1.01 ng m⁻³). Melpitz shows similar retene concentrations to Frýdlant (1.16 ng m⁻³) at lower total PAH levels (14.12 ng m⁻³). The high relative and absolute levels of retene and levoglucosan in Košetice are a good indicator of a high contribution of wood burning. The results are also consistent with those of Horník et al. (2024), who reported high levels of levoglucosan in Košetice and Frýdlant. Overall, the results indicate a strong influence of wood burning for domestic heating during winter in the Košetice area. During particularly cold periods, residents appear to supplement wood with coal, leading to a greater local impact of coal emissions on air quality. The higher coal contributions observed in Melpitz seem to be mainly driven by long-distance transport, whereas in Frýdlant, additional contributions from the nearby power plant are evident.

4.1 Effects of COVID-19 containment measures

With containment measures still in place during the campaign, the daily lives of many of the region's citizens were disrupted. Different patterns of mobilisation, the closure of businesses, and changes in leisure habits are all factors that can affect air quality. The emission inventories used in this study do not take into account exceptional events affecting emissions, such as the COVID-19 restrictions. Several studies have looked at the impact of these restrictions on air quality. Most of them focus on the year 2020, when the pandemic peaked.

Gkatzelis et al. (2021) reviewed over 200 papers to assess the impact of lockdowns on air quality around the world. They found significant reductions in NO₂ and CO levels, small reductions in PM_2.5, and increases in O₃ concentrations. The effects varied by season and region, and the study highlighted the need for future research to include meteorological corrections for accurate results. Only about a third of the studies reviewed included methods for meteorological correction or normalisation.

Matthias et al. (2021) conducted a modelling study for central Europe, estimating emission reductions from January to June 2020. For secondary pollutants, they found that meteorological effects outweighed the effects of emission reductions from restrictions. Putaud et al. (2023) compared measurements at 28 sites across Europe for spring 2020 with CAMS ensemble forecasts and found a slight decrease in PM_2.5 and PM₁₀ during the lockdown and a strong increase after the measures were lifted. The study corrects the bias occurring between modelled and measured values through a time-dependent normalisation of the CAMS forecasts to the observations estimated from 2019 data. They concluded that the increased ozone levels due to reduced NO_x led to altered oxidation capacities and therefore more SOA formation. The study also analysed data collected in Melpitz and Košetice before, during, and after the lockdown in March 2020. In Melpitz, slightly higher PM_2.5 concentrations than expected by CAMS were detected during the lockdown. In May 2020, after the lockdown, they were even twice as high as the modelled concentrations. In Košetice, the values before and during the lockdown were slightly below the expected values, while the concentrations afterwards were 30 % higher. Forster et al. (2020) calculated emission trends based on Google mobility data for six sectors (land transport, residential, energy, industry, public, and aviation) per country. These data show that in March 2021, BC emissions from the residential sector in Germany and the Czech Republic were increased by approximately 10 %, while BC emissions in all sectors combined were decreased by about 20 % compared to a baseline scenario. Mbengue et al. (2023) conducted an extended study analysing the effects of the COVID-19 lockdowns in Košetice using normalisation techniques to account for meteorological effects. They found that during the winter of the second lockdown (December 2020–February 2021), dispersion normalised concentrations of EC were reduced by 28 %, while OC and SOC concentrations increased by 19 % and 51 %, respectively. They concluded that this was due to a greater influence of emissions from local domestic activities. Considering that our study sites are background stations with a low traffic influence and a high contribution of domestic heating emissions, locally increased emissions due to the COVID-19 mitigation measures seem plausible, leading to higher PM_2.5 and probably SOA concentrations than without these measures. These changes are not included in the emissions in the model and may be another source of underestimation in the model.

4.2 Anthropogenic secondary organic aerosol

Bergström et al. (2012) found an underestimation of winter organic aerosol in a modelling study focusing on several years in Europe. Their conclusion was that emissions from wood combustion are underrepresented in current emission inventories. Previous source apportionment studies have shown that residential heating is a significant contributor to SOA formation. Lhotka et al. (2025) identified a relationship between primary organic aerosol (POA) and oxidised organic aerosol (OOA) source factors associated with residential heating. The high contribution of highly oxidised OA in winter can be attributed to the local influence of biomass burning. In contrast, in Melpitz, coal combustion plays a more prominent role in oxidised OA formation, indicating the impact of long-range transport (Atabakhsh et al., 2023). An intensive tagging study by Bartík et al. (2024) utilised the PSAT module in CAMx, supplementing the CAMS emission inventory with more detailed residential emission data for the Czech Republic and additional intermediate-volatility organic compound (IVOC) emissions from wood combustion. Their findings indicate that VOC and IVOC emissions from the “other combustion” sector represent the largest source of SOA in central Europe during winter, contributing up to 0.4 µg m⁻³. In order to investigate whether a potential underestimation of SOA precursors from domestic heating contributed to the lower-than-expected concentrations of OM in our model, we carried out a sensitivity study.

4.2.1 Sensitivity study

The parameterisation of SOA is influenced by two key variables: the precursor gases emitted and the rate at which SOA is formed from these precursors. Previous studies suggest that phenol is a significant component of emissions from incomplete combustion processes like wood burning. Phenol is one of the key gaseous precursors responsible for the formation of SOA during biomass burning activities (Hatch et al., 2015; Bruns et al., 2016). Liu et al. (2024) highlight the critical role of nighttime NO₃ oxidation of anthropogenic VOCs from biomass combustion, a process they find is inadequately represented in current atmospheric models. Their results show that increasing both phenol emissions and the associated SOA yield leads to a 2-fold increase in SOA production via NO₃ oxidation across Europe during winter. In our model, phenol is included in the lumped species cresol and other aromatics (CSL) (see Reaction 1). NMVOC emissions are delivered by the UBA and CAMS emission inventories (Schneider et al., 2016; Kuenen et al., 2022). The NMVOC emission flux is split into the different relevant model species based on emission profiles created by Theloke and Friedrich (2007) for 306 individual species, including phenol. These profiles are based on an NMVOC source database from 1990 (Olivier et al., 1996) and do not include phenol emissions from domestic heating. Therefore, it is reasonable to assume that the VOCs emitted from domestic heating are not fully captured by the model. Natural fire emissions provided by GFAS do not include CSL but rather toluene (TOL) and xylene (XYL) emissions. Given the low impact of natural fires in winter in Europe, it can be assumed that they do not contribute much to the formation of secondary particles.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{CSL}+\mathrm{OH}& \to {\mathit{\alpha }}_{\mathrm{1}}\mathrm{CVARO}\mathrm{1}+{\mathit{\alpha }}_{\mathrm{2}}\mathrm{CVARO}\mathrm{2}\\ \mathrm{CSL}+{\mathrm{NO}}_{\mathrm{3}}& \to {\mathit{\alpha }}_{\mathrm{1}}\mathrm{CVARO}\mathrm{1}+{\mathit{\alpha }}_{\mathrm{2}}\mathrm{CVARO}\mathrm{2}\end{array}\end{array}$$

The SORGAM module (see Sect. 2.3) estimates SOA formation from aromatic precursors using data from Odum et al. (1997), who conducted smog chamber experiments to quantify SOA production from gasoline vapour. This method is therefore primarily tailored to traffic emissions. However, for aromatic precursors emitted by domestic heating, an increase in SOA yield aligned with phenol SOA formation rates is more suitable. Due to the limited availability of chamber studies on phenol gas-phase SOA formation, we derived a new yield estimate based on four OH oxidation measurements from Yee et al. (2013). Given the importance of nocturnal oxidation, we also applied these modifications to the NO₃ reaction. A non-linear least squares fit for the α_i values was performed with fixed-K_om,i coefficients ($K_{\mathrm{om},\mathrm{1}}=\mathrm{0.2899}$, $K_{\mathrm{om},\mathrm{2}}=\mathrm{0.0103}$). As α₂ yielded negative values, we decided to keep α₂ fixed and performed the fit solely for α₁ (see Fig. A6 in the Appendix). These adjustments result in an approximately 3-fold increase in SOA yield.

Further, an adjustment was made to the input emission. To obtain a good representation of phenol emissions from domestic heating processes, we decided to scale the emissions to the CO emissions of the emission sector “other combustion”. Following wood combustion chamber studies from Bruns et al. (2016), on average all NMVOC emissions make up 0.22 times the CO emissions. According to Schauer et al. (2001), phenol and substituted organic compounds are approximately 10 % of the overall NMVOC emissions from wood combustion. Accordingly, we set our “other combustion” sector CSL emissions to 0.022 times the sector's CO emissions.

We simulated three sensitivity runs to compare these adjustments: first with the adjusted SOA yield alone (S1), second with the new CSL emissions alone (S2), and third with both combined (S3). Table 4 gives an overview of the coefficients and emissions used in the different sensitivity runs and the original base run. All sensitivity runs were performed for our middle domain, TRACED0 (see Fig. 1), as it provides the best trade-off between spatial resolution and area coverage. The three sensitivity runs were not nested but used the same initial and boundary conditions as the base run.

Table 4 Overview of the sensitivity simulations. Shown are changes to α₁ to adjust the SOA yield parameter for aromatic precursors and scaling of CSL emissions based on CO emissions from GNFR C to account for phenol contributions.

4.2.2 Sensitivity study results

The changes in CSL emission flux and the corresponding mean OM concentration across the three sensitivity runs are presented in Table 5. The values for each station represent the model result from the 4×4 km grid cell in which the station is located. In S1, the increased SOA yield for aromatic precursors has the most pronounced effect in urban areas, as it influences emissions from all source sectors, including industry and transport. The adjusted SOA yield applies to both daytime OH oxidation and nighttime NO₃ oxidation. Among the stations, Melpitz shows the highest relative increase (39 %) in mean OM concentrations, reaching 1.86 µg m⁻³, due to high aromatic precursor levels.

Table 5 Changes in CSL emission flux and mean organic matter in PM_2.5 for all sensitivity runs compared to the base run.

In S2, the emissions of aromatics from domestic heating are introduced as CSL emissions by the sector “other combustion”. Although total CSL emissions across the domain remain constant, their spatial distribution shifts: emissions decrease in Melpitz but increase significantly in Frýdlant and Košetice. The domain-wide mean OM concentration shows an overall modest increase of 4 %, with the largest increases observed in the central Czech Republic and southern Poland, where domestic heating sources are abundant. Interestingly, despite a reduction in CSL emissions in Melpitz compared to the base run, OM concentrations in Melpitz increase similarly to those in Frýdlant (+16 % in Melpitz and +15 % in Frýdlant). This is attributed to increased CSL emissions in the surrounding areas and the transport of SOA and its precursors to the site. These findings align with previous studies: Poulain et al. (2011) linked winter OM in Melpitz to transported particles, while Spindler et al. (2012) reported that SOA concentrations peaked in winter air masses arriving from the east, highlighting the role of anthropogenic-precursor-driven SOA formation during long-range transport. This is also consistent with the conclusions of Atabakhsh et al. (2023).

In the combined S3 run, average OM concentrations in Melpitz increase by 58 % to 2.11 µg m⁻³, representing the highest relative impact among all stations. This increase can be attributed to enhanced SOA transport and formation from aromatic precursors. Frýdlant shows the largest absolute OM increase, with an average increment of 0.95 µg m⁻³, reaching 2.76 µg m⁻³ (see Table 5). Figure 9 compares diurnal OM cycles from the base and sensitivity runs with measurements. In Frýdlant and Košetice, the combined adjustments in S3 produce greater impacts on OM concentrations compared to the individual sensitivity runs. This leads to better agreement with the measurements in Košetice but results in overestimation compared to the AMS/ACSM data in Frýdlant. As discussed previously, the discrepancies between the AMS/ACSM and Sunset Laboratory measurements cannot be fully resolved in this work, and both datasets must therefore be regarded as valid. Taking into account the measurement uncertainties, the fact that the simulated OM concentrations in Frýdlant now lie between the two measurements supports the plausibility of the modelled increase. Evaluating both datasets in combination provides a more comprehensive and balanced assessment of actual OM levels. The AMS/ACSM is better suited to capturing diurnal patterns due to its higher time resolution. In Frýdlant, the model simulates a clear morning peak in OM concentrations that is absent in the AMS/ACSM data. This discrepancy suggests that the model may be overestimating the contribution from local or near-field sources, while underestimating the influence of long-range transport.

In Melpitz, however, the difference between S1 and S3 is less distinct, suggesting that higher baseline precursor concentrations already contribute significantly to SOA formation at this site. The correlation of the modelled SOA with AMS/ACSM data improves in Košetice, with the correlation coefficient increasing from 0.18 to 0.29, while Melpitz and Frýdlant show no significant improvement. Although the model now better reflects SOA contributions in Košetice, overall OM concentrations remain underestimated.

Figure 9 Hourly diurnal cycle of organic matter in PM_2.5 over the entire measurement period for all sensitivity runs compared to the base run. Bar graphs show averaged filter data based on 12 h sampling intervals, with whiskers representing a measurement uncertainty of ±12 % (Karanasiou et al., 2020). AMS/ACSM data are shown with a measurement uncertainty of ±25 % (Canagaratna et al., 2007).

Figure 10 Difference between the averaged modelled organic matter in PM_2.5 in the S3 and base simulation.

The spatial pattern of the increase in OM concentrations in S3 compared to the base run is shown in Fig. 10. The strongest increases occur in the Czech Republic and southern Poland. In particular, the city of Prague and its outskirts show a high increase, which is consistent with the small residential houses in the suburbs and surrounding areas still using coal and wood combustion for heating (Schwarz et al., 2008; Makeš et al., 2021). Domestic heating emissions likely increased during the COVID-19 measures as more people stayed at home, contributing to higher PM_2.5 and possibly SOA concentrations (Mbengue et al., 2023). Additionally, the campaign coincided with an unusually cold period in early February, which likely increased heating activity further. This effect was not fully captured by the model.

The simulation of the OM peak in early February, after the snow event, shows no noticeable improvement with the sensitivity runs, likely linked to the overall underestimation of emissions (see Fig. A8 in the Appendix). Measurements reveal a distinct OM concentration peak on 3 March, particularly in Košetice and Frýdlant. The S3 run captures the peak reasonably well in Frýdlant but still underestimates OM in Košetice. HYSPLIT backward trajectories for Frýdlant on this date indicate a significant influence from air masses passing over the Czech Republic and Košetice (see Fig. A7 in the Appendix). In Košetice, the wind direction shifts from east to west, as simulated by the model, while in Frýdlant, the wind direction remains steady, allowing precursor accumulation and increased SOA formation.

Overall, the sensitivity studies showed that scaling AVOC emissions from wood combustion to residential heating emissions improves the spatial distribution of SOA in the study area. Long-range transport of precursors and SOA is captured as well as the local influence on OM concentrations.

5 Conclusion, limitations, and future directions

The study investigates the sources of primary and secondary anthropogenic organic aerosol in central Europe during winter. The chemical transport model COSMO–MUSCAT was used to analyse concentrations of particulate matter, in particular particles originating from combustion processes. The model results were compared with measurements made in Germany and the Czech Republic in terms of overall PM_2.5 concentration and concentrations of individual species. The model underestimated the total PM_2.5, especially during high concentration peaks. A pronounced underestimation occurred in early February, likely due to the prevailing meteorological conditions combined with changed heating behaviour. During this period, all tracers were underestimated, whereas after early February, the model accurately captured the behaviour of most tracers. However, the underestimation of PM_2.5 during concentration peaks remained evident. The discrepancies in modelled PM_2.5 concentrations do not appear to be due to deviations in EC, sulfate, or nitrate levels but rather the underestimation of OM. Although the present study reproduced diurnal OM profiles well at two monitoring sites, measurements in Košetice are underestimated, partly due to an inadequate representation of SOA formation from residential heating (wood combustion), a major source of anthropogenic VOCs. These AVOCs contribute considerably to the formation of SOA, and it is likely that their insufficient representation in our model contributes to the overall underestimation of OM during winter. The effect is most pronounced in the central Czech Republic, where the basin-like topography allows air masses to linger, promoting the accumulation of emissions and extended SOA formation. We found a higher contribution of domestic heating in the eastern part of our study region, which is accompanied by high concentrations of OM, especially at the station in Košetice. Sensitivity tests with adjustment for SOA yields and AVOC emissions showed an average increase in OM concentrations of over 40 % at the measurement sites. It is likely that the model underestimates SOA precursor emissions from domestic heating sources and additional sources that are missing or unaccounted for in the underlying emission inventory. A more detailed inventory, as used in Bartík et al. (2024) for the Czech Republic, reveals a redistribution of total primary PM_2.5 residential combustion emissions from urban to rural areas compared to the inventory used in this study. In addition, the model may underestimate the contribution of SOA precursors other than phenol. Implementing more detailed and up-to-date emission inventories that provide information on the types of fuels used, their spatial distribution, and their temporal profiles offers strong potential to enhance the model's overall performance with respect to OM. Our findings highlight that regional domestic heating emissions contribute significantly to overall air pollution in the study area. Addressing these emissions is complex, as they are hard to quantify, and regulations for private households are more challenging to implement. Consequently, obtaining more detailed information on these sources is vital for developing targeted and feasible measures. Besides updated time profiles representing seasonal, weekly, and daily patterns of emissions, changing heating behaviour due to extreme meteorological conditions could be taken into account by implementing a temperature dependence of emissions. The heating degree day (HDD) approach, introduced by Guevara et al. (2021), considers the influence of outdoor temperature on heating activity and its associated emissions. A recent study by Guion et al. (2024) enhanced this method by incorporating country-specific and species-specific parameters, demonstrating improved temporal correlations and more accurate detection of PM emission threshold exceedances compared to simulations using fixed parameters or monthly temporal factors. Implementing this approach in COSMO–MUSCAT could enhance the accuracy of our model results during winter. In addition, heating emissions may not only be underestimated in quantity, but the contribution of different fuel types to the domestic heating sector may also vary with temperature, as additional coal burning in households may occur during colder periods.

Comparing the non-reactive tagging approach in COSMO–MUSCAT to measurement-based, receptor-oriented source apportionments can further help evaluate its capability and identify areas for improvement. This comparison can provide valuable insights into the performance of the model and guide future refinements. By addressing these gaps and incorporating the necessary updates, such as updated emission inventories, improved SOA yields, and model evaluation through comparison with measurement data, the model could provide a more comprehensive and accurate representation of SOA formation processes, enabling a better understanding for air quality management.

Appendix A

Figure A1 Time profile applied to the GNFR C emission factor following Kuenen et al. (2014).

Table A1 Splitting factors applied in this study to disaggregate CAMS PM_2.5 and PM₁₀ bulk emissions into individual subgroups within the GNFR C emission sector (Kuenen et al., 2022).

Table A2 Splitting factors applied in this study to disaggregate CAMS NMVOC bulk emissions into individual subgroups within the GNFR C emission sector (Kuenen et al., 2022).

Table A3 Overview of all the tagged combinations of variables, source regions, and source sectors.

Figure A2 Left: time series for organic matter mass concentration for the three sites. Comparison of measured AMS/ACSM data with the modelled primary OM concentration, split into the share with origin outside the N0 domain, origin inside the study domains, and secondary organic aerosol. Right: normalised POA to total OM ratio.

Figure A3 Time series for elemental carbon concentration for the three sites. Comparison of Aethalometer and Sunset Laboratory filter data and modelled data.

Figure A4 Time series for nitrate mass concentration for the three sites. Comparison of measured and modelled data.

Figure A5 Time series for sulfate mass concentration for the three sites. Comparison of measured and modelled data.

Figure A6 Non-linear least squares fit for two-product aromatic class.

Figure A7 Backward trajectory ending on 3 March 2021 at 12:00, 15:00, and 18:00 UTC in Frýdlant, created with the NOAA HYSPLIT trajectory model (Rolph et al., 2017; Stein et al., 2015).

Figure A8 Time series for organic matter mass concentration for all sensitivity runs compared to the base run and measurements.

Table A4 Relative contributions of different source sectors and source regions to the total of ECfine, ECcoarse, OMfine, and OMcoarse and absolute mean over all sectors. Contributions from outside the N0 domain are not included.