Measurement report: Source attribution and estimation of black carbon levels in an urban hotspot of the central Po Valley – an integrated approach combining high-resolution dispersion modelling and micro-aethalometers

Understanding black carbon (BC) levels and its sources in urban environments is of paramount importance due to the far-reaching health, climate, and air quality implications. While several recent studies have assessed BC concentrations at specific fixed urban locations, there is a notable lack of knowledge in the existing literature on spatially resolved data alongside source estimation methods. This study aims to fill this gap by conducting a comprehensive investigation of BC levels and sources in Modena (Po Valley, Italy), which serves as a representative example of a medium-sized urban area in Europe. Using a combination of multi-wavelength micro-aethalometer measurements and a hybrid Eulerian–Lagrangian modelling system, we studied two consecutive winter seasons (February–March 2020 and December 2020–January 2021). Leveraging the multi-wavelength absorption analyser (MWAA) model, we differentiate sources (fossil fuel combustion, FF, and biomass burning, BB) and components (BC vs. brown carbon, BrC) from micro-aethalometer measurements. The analysis reveals consistent, minimal diurnal variability in BrC absorption, in contrast to FF-related sources that exhibit distinctive diurnal peaks during rush hours, while BB sources show less diurnal variation. The city itself contributes significantly to BC concentrations (52 ± 16 %), with BB and FF playing a prominent role (35 ± 15 % and 9 ± 4 %, respectively). Long-distance transport also influences BC concentrations, especially in the case of BB and FF emissions, with 28 ± 1 % and 15 ± 2 %, respectively. When analysing the traffic-related concentrations, Euro 4 diesel passenger cars considerably contribute to the exhaust emissions. These results provide valuable insights for policy makers and urban planners to manage BC levels in medium-sized urban areas, taking into account local and long-distance sources.

1 Introduction

Both black carbon (BC) and elemental carbon (EC) are carbonaceous particles associated with particulate matter (PM) that result from incomplete combustion processes and have light-absorbing properties. Although they are commonly mistaken for synonyms, these two terms should not be treated as completely interchangeable. BC denotes materials primarily composed of carbon with high light-absorption potential, while EC specifically denotes pure carbon without bonds to other elements (Petzold et al., 2013). While these formal definitions strive for rigour, their practical application can be challenging due to the inherent variability in carbonaceous matter. Unlike pure substances, real-world BC exists as a complex mixture of various compounds with different properties. For this reason, a common approach to defining BC is based on properties such as light absorption, solubility, thermal stability, morphology, or microstructure (Petzold et al., 2013).

BC originates from a variety of sources, including combustion engines (especially diesel), residential burning of wood and coal, and field burning of agricultural residue, as well as forest and vegetation fires. The detrimental consequences associated with BC have become a key element of scientific investigation and public concern, necessitating a comprehensive understanding of its sources, behaviour, and potential mitigation strategies. The adverse effects of BC on human health are well-documented in the literature. Inhalation of BC particles has been linked to various respiratory and cardiovascular conditions, including asthma, chronic bronchitis, reduced lung function, and increased risk of heart attacks and strokes (Song et al., 2022; Janssen et al., 2012; Grahame et al., 2014; Rohr and McDonald, 2016; Janssen et al., 2011). Moreover, long-term exposure to BC has been associated with an elevated incidence of lung cancer (Chang et al., 2022). These health risks are particularly worrisome considering the widespread prevalence of BC in urban areas (Reche et al., 2011; Ali et al., 2021) and its ability to penetrate deep into the respiratory system (Saputra et al., 2014; Tang et al., 2020), posing a significant threat to vulnerable populations, such as children, the elderly, and individuals with pre-existing respiratory disorders.

BC not only poses adverse health effects but also exerts significant impacts on the Earth's radiation balance through various mechanisms (Chung and Seinfeld, 2005; Wang, 2004; Roberts and Jones, 2004). These include direct and indirect effects, altering the amount, distribution, and properties of solar energy (Team et al., 2023). The direct effects encompass the absorption and scattering of sunlight by BC, leading to localised atmospheric warming and influencing the distribution of solar radiation within the atmosphere (Menon et al., 2002; Ramanathan et al., 2001; Ramanathan and Carmichael, 2008). Estimates of BC's direct radiative forcing range from +0.1 to +1.0 W m⁻², exhibiting variations depending on the methodology employed, such as bottom-up emission inventories and atmospheric transport model approaches or observation-based estimates (Wang et al., 2016; Bond et al., 2013). In addition to the direct effects, BC particles have indirect impacts on the Earth's radiation balance by interacting with clouds. BC can serve as cloud condensation nuclei (CCN) or ice nuclei (IN), influencing the formation of cloud droplets and ice particles (Hendricks et al., 2011; Koch et al., 2011). The presence of BC can modify cloud microphysical properties, including droplet size, concentration, cloud lifetime, and precipitation processes. These indirect effects have been estimated to range from −0.13 to −0.11 W m⁻², contributing to negative radiative forcing (Cherian et al., 2017; Koch et al., 2011). Considering both the direct and indirect effects, the combined impact of BC on the Earth's radiation balance results in positive radiative forcing on the climate system, exacerbating global warming and climate change.

In the literature, several approaches have been employed to measure BC, EC, or more generally, light-absorbing aerosols. These approaches can be categorised into four main methods: thermo-optical determination, photothermal interferometry, photoacoustic spectroscopy, and optical determination. Thermo-optical determination involves the measurement of light-absorbing aerosols by analysing their thermal properties. The technique typically relies on heating the sample in a controlled environment and measuring the resulting changes in optical properties, such as light absorption and scattering, as a function of temperature. By observing these changes, one can estimate the amount and type of light-absorbing aerosols present in the sample (Chow et al., 2007; Bauer et al., 2009; Brown et al., 2019). Field experiments adopting this determination method are, for example, Merico et al. (2019), Liakakou et al. (2020), and Bigi et al. (2017).

Photo-thermal interferometry is a technique that combines laser-induced heating with interferometric measurements to quantify light-absorbing aerosols (Lee and Moosmüller, 2020; Visser et al., 2023; Li et al., 2016). By directing a laser beam onto the aerosol sample, the absorbed light generates local temperature gradients, causing changes in the refractive index and thus altering the phase of a probe laser beam passing through the sample. These phase changes are then measured interferometrically to determine the concentration of light-absorbing aerosols. On the other hand, photo-acoustic methodologies involve the use of laser-induced acoustic waves to detect and quantify light-absorbing aerosols. In this technique, a pulsed laser is used to heat the sample, causing the aerosols to absorb the light and generate acoustic waves. By measuring the amplitude of these acoustic waves, the concentration of light-absorbing aerosols can be determined (Fischer and Smith, 2018; Petzold and Niessner, 1996; Guo et al., 2014). The deployment of photo-thermal and photo-acoustic instruments in field settings have been limited due to operational challenges, which has hindered the widespread adoption of these techniques. However, there has been a recent resurgence of interest in these methods with the development of a new generation of photometers and interferometers (Drinovec et al., 2022; Visser et al., 2020, 2023).

In addition to the aforementioned techniques, optical determination has gained popularity in recent years for measuring the light absorption properties of BC aerosols due to its cost-effectiveness, comparability with chemical analysis, and ability to provide high temporal resolution results for real-time monitoring studies (Kaskaoutis et al., 2021a, b; Yus-Díez et al., 2021; Bernardoni et al., 2021). In this approach, sample air is passed through a filter tape where aerosol particles are collected. Optical filter photometers measure light transmission, or a combination of reflection and transmission, through the sample-loaded filter and calculate the attenuation coefficient from the rate of change in attenuation over time. However, the attenuation coefficient can differ significantly from the true aerosol absorption coefficient due to two main artefacts. The first is the enhancement of the optical path, and hence the enhancement of light absorption of the deposited particles, due to the multiple scattering of the light beam at the filter fibres and between particles and fibres. The second is the loss of instrument sensitivity with increasing particle loading. To overcome these limitations, several empirical corrections have been proposed in the literature. For instance, the C_ref factor is used to correct for the multiple scattering effect, while the f(ATN) function is applied to compensate the loading effect. Typically, the C_ref factor is assumed a priori, but it can also be determined experimentally through multi-instrument colocation. Conversely, the loading correction function f(ATN) is increasingly estimated online using dual-spot technology (Drinovec et al., 2015), or it can be estimated offline using dedicated algorithms (Weingartner et al., 2003; Virkkula et al., 2007; Park et al., 2010). The specific mass absorption efficiency (MAE; also known as the mass absorption cross-section, MAC) can then be used to convert the aerosol light absorption coefficient into the light-absorbing carbon mass concentration. Although this relationship is simple, the value of the MAE varies considerably in time and space depending on emission sources, transport phenomena, combustion conditions, particle ageing, and mixing state (Chan et al., 2011; Mbengue et al., 2021), making the conversion process a source of potential uncertainty (Petzold et al., 2013). When considering BC material, the term equivalent BC (eBC) is commonly used to report the mass concentration indirectly determined by light absorption techniques such as filter-based absorption photometers (Petzold et al., 2013). Therefore, the term eBC will be used throughout this study to refer to those compounds calculated using filter-based absorption photometers.

BC exhibits distinctive light absorption characteristics, particularly in the infrared range, facilitating its identification and quantification. Widely used instruments exploiting this technique are the aethalometer, such as the AEs models (e.g. AE31, AE33, AE36, and AE43; Magee Scientific Co.) or the portable micro-aethalometer MA series (e.g. MA200, MA300, and MA350; Aethlabs), which are filter-based optical instruments operating at different wavelengths. This later feature, coupled with the usage of source-specific absorption Ångström exponent, can be used for the source apportionment of black carbon (Zotter et al., 2017; Sandradewi et al., 2008; Massabò et al., 2015; Bernardoni et al., 2017b).

Source apportionment of BC is essential to gain insight into the specific pollution sources contributing to atmospheric concentrations. This analysis enables the development of targeted mitigation strategies, in line with the recommendation proposed by the EU legislation (European Council, 2008), which highlights the importance of identifying and addressing specific emission source sectors to effectively reduce air pollution. Furthermore, the World Health Organization (WHO), in its recently updated global air quality guidelines, has recommended the systematic measurement of BC (or EC), the establishment of BC inventories, and, where appropriate, the implementation of BC reduction measures (WHO, 2021). In 2022, the European Commission proposed revisions to the Ambient Air Quality Directives to bring European Union air quality standards more closely in line with WHO recommendations (European Council, 2022). These revisions highlight the importance of monitoring emerging pollutants such as BC to support scientific understanding of their effects on health and the environment, in line with the WHO's guidance. This research aims to provide information on the levels of BC and to identify the main sources of BC in the urban environment of Modena, a city of about 200 000 inhabitants located in the middle of the Po Valley, a densely populated and industrialised region in northern Italy known for its high levels of air pollution (Bigi and Ghermandi, 2016; Bigi et al., 2012; Lonati et al., 2017; Perrino et al., 2014; Veratti et al., 2023; Thunis et al., 2021).

Previous studies have investigated the sources and transport of particulate matter in the Po Valley using both chemical transport models (CTMs) and aerosol composition, such as, for example, Scotto et al. (2021), Paglione et al. (2020), Bernardoni et al. (2011), Pepe et al. (2019), Belis et al. (2019), and Bernardoni et al. (2017a), but there is still a knowledge gap regarding the contribution of different sources to BC concentrations. For instance, a study by Mousavi et al. (2019) combined the usage of an aethalometer, thermo-optical instruments, and ¹⁴C analysis to apportion BC between fossil fuel and biomass burning for three sites in the urban area of Milan and its surroundings. Similar studies in the Po Valley have been conducted employing a variety of instruments such as aethalometers, multi-angle absorption photometers (MAAP), polar-photometers, and independent measurement of levoglucosan (Bernardoni et al., 2021; Massabò et al., 2015; Bernardoni et al., 2017b; Gilardoni et al., 2020). In addition, concerns about health effects have led recently to a surge of research in urban areas, highlighting eBC concentrations as a fundamental tracer of pollution in cities (Segersson et al., 2017; Pani et al., 2020). As a result, several studies have investigated the variability in eBC mass concentrations and the associated source apportionment at multiple urban locations. Examples of this research in Europe are Savadkoohi et al. (2023), Liakakou et al. (2020), Grange et al. (2020), Helin et al. (2018), Minderytė et al. (2022), Titos et al. (2017), and Singh et al. (2018). Although all these studies utilised accurate methods for eBC determination and source apportionment, their results were limited to specific monitoring sites, including rural or urban background locations. As a consequence, they were unable to provide continuous spatial information across the entire city, which is important information for consistent epidemiological analysis.

In addition to the aforementioned studies, other researchers have conducted modelling studies at the European level to assess the performance of air quality models in reproducing carbonaceous aerosol concentrations and their absorption properties, such as Mircea et al. (2019) and Curci et al. (2019). Other work aimed to estimate BC concentrations at the street level in small neighbourhoods in Paris, Brussels, or Maribor using a state-of-art multi-scale modelling system, Street-in-Grid (Lugon et al., 2021), or a more simplified approach (Brasseur et al., 2015; Ježek et al., 2018). However, none of these studies focused on assessing total BC concentrations over an entire urban area at high resolution or attempted to identify the contribution of different sources at monitoring sites.

In this study, we conducted a two-sampling-period campaign during wintertime (February–March 2020 and between December 2020 and January 2021) in Modena, deploying two MA200 multi-wavelength micro-aethalometers at distinct urban locations and at traffic and background sites. The primary objectives were to determine total BC concentrations at both sites and to investigate the potential sources contributing to the observed concentrations. To support the identification of local and regional sources and to provide spatial information about BC concentrations, we employed a hybrid Lagrangian–Eulerian modelling system composed of GRAMM-GRAL (the Graz Mesoscale Model – Graz Lagrangian Model; Oettl, 2015b, a, c) and the NINFA (Network dell'Italia del Nord per previsioni di smog Fotochimico e Aerosol; Vitali et al., 2023) modelling suite, which is built upon the CHIMERE chemical transport model (Mailler et al., 2017; Menut et al., 2021). By combining observed and modelled data, we aimed to improve our understanding of the main emission sectors contributing to total BC levels for the development of effective mitigation strategies and aimed to provide valuable insights for addressing environmental and health concerns related to BC.

The paper is structured as follows. Section 2 provides an overview of the measurement campaign and the models that make up the hybrid system. Further details on emission estimation are given in Sect. 2.5. Section 3 presents the statistical metrics used to validate the model performance against the observational data. Section 4 presents the results of the measurement campaign and modelling activities. Finally, Sect. 5 presents some concluding remarks.

2 Materials and methods

2.1 Sampling sites and measurement description

The measurement campaign was carried out in Modena, a medium-sized city in the Emilia–Romagna region of northern Italy. Modena is located about 50 km south of the Po River in the Po Valley, an area known for its dense population and industrialisation. The city, like the whole valley, has a continental climate with marked seasonal variations. In particular, during the winter months, the Po Valley has specific meteorological characteristics, such as stable atmospheric conditions and frequent low wind speeds, which contribute to the occurrence of thermal inversions. These atmospheric phenomena can hinder the dispersion of air pollutants in the region, as reported by many studies, such as Bigi et al. (2012), Pernigotti et al. (2012a), Pernigotti et al. (2012b), and Ghermandi et al. (2017). In addition, the Po Valley is geographically bounded by the Apennines to the south and the Alps to the north, which act as natural barriers that further limit the movement of pollutants. As a result, the combination of high traffic volumes, industrial activities, and the large number of densely populated urban areas in the valley can lead to the accumulation of pollutants.

In this study, aerosol measurements were performed using two MA200 micro-aethalometers (microAeth^® MA Series, Aethlabs, USA) equipped with PTFE filter tapes. These instruments quantify the absorption of light-absorbing carbonaceous aerosols at five different wavelengths (375, 470, 528, 625, and 880 nm), allowing for the estimation of eBC concentrations by employing wavelength-specific mass absorption efficiencies (MAE).

Sampling was performed at two air quality monitoring stations in the city, utilising the existing glassware manifold inlet lines used for reactive-gas monitors. The inlet temperature was maintained at 30 ± 2 °C, and no size cut-off was applied to the aerosol particles entering the devices. One station is located in a public park on the west side of the city, representing typical background conditions where pollution levels are influenced by a mixture of contributions rather than a single source type, while the other station is located in an urban traffic area characterised by high vehicular activity, serving as a representative location for areas of dense traffic and associated air pollution. The locations of the two stations are shown in Fig. 1.

Figure 1 Spatial representation of the NINFA and GRAMM domains on the left (from Esri, USGS, NOAA) and the GRAL investigation domain on the right. The GRAL domain includes the positions of two urban air quality stations denoted by light-blue dots, specifically one located at the urban traffic (traffic) site and another at the urban background (background) site. Additionally, four meteorological stations used for meteorological validation and for selecting the mesoscale meteorological situation are indicated by yellow stars within the GRAL domain.

Throughout the study, the instruments operated at a flow rate of 100 mL min⁻¹ in dual-spot sampling mode to compensate for filter loading effects (Drinovec et al., 2015). Although several papers in the literature have emphasised the need for a multi-instrument approach to determine the multi-scattering correction factor (C_ref) for aethalometer filters (Bernardoni et al., 2021; Yus-Díez et al., 2021; Ferrero et al., 2021), with photothermal interferometry having recently been identified as a suitable technique to provide a reference measurement of this parameter (Drinovec et al., 2022), this study lacked multi-instrument co-location. As a result, we opted for a constant C_ref of 1.3 as suggested by the manufacturer to mimic the response of the AE33 aethalometer (Aethlabs, 2024). In addition, to convert the aerosol light absorption coefficient to the equivalent mass concentration, we relied on the wavelength-specific MAE values provided in the MA200 reference manual and reported in Table S1 in the Supplement. In support of these decisions, a recent instrument intercomparison performed at an urban background site in Athens (Stavroulas et al., 2022) showed limited differences in terms of eBC concentrations between an MAAP and the two MA200 units used in this study, set to a C_ref of 1.3 and default MAE values (linear slope of 1.00 in winter and 1.07 in summer, with an r² of 0.92 for both the seasons). From the same intercomparison campaign, two MA200 units were also compared with an AE33 aethalometer, showing strong agreement during winter (linear slope between 0.91 and 0.97 and r² of 0.97 for both devices). Furthermore, recent studies corroborate the results of the Athens intercomparison, showing consistent findings when comparing the AE33 and the MA200. For example, Blanco-Donado et al. (2022) observed a linear regression slope of 0.97 and a r² of 0.93 during a 3 d intercomparison campaign in a suburban area of Barranquilla (Colombia). Similarly, Khan et al. (2024) reported a linear regression slope of 0.986 and a r² of 0.97 for a 14 h intercomparison conducted in an urban background area in Vilnius (Lithuania).

Table S1 presents a comparison of the MAE values used in this study with those measured across various European locations. Although the values used for Modena are generally higher than those reported in the literature, they are consistent with previous observations conducted in the Po Valley, such as Gilardoni et al. (2020) and Mousavi et al. (2019). Furthermore, the MAE value of 880 nm used for this study is in the range of values reported for other European urban areas, such as in Athens (Savadkoohi et al., 2024), Zurich, and Bern (Grange et al., 2020), or in the range observed for urban traffic and background sites in Leipzig and Prague for MAE values of 637 nm (Savadkoohi et al., 2024).

Measurements at the urban traffic site were made during two periods: from 4 February to 7 March 2020 and from 26 December 2020 to 21 January 2021 with a time resolution of 1 min, while the background site was monitored from 4 February to 7 March 2020 with a time resolution of 1 min and from 26 December 2020 to 7 January 2021 with a time resolution of 5 min. To account for occasional low absorption readings at the background site, the 1 min data at the urban background location were aggregated to 5 min using the dual-spot compensation algorithm implemented in the R language (Bigi et al., 2023), following the formulation proposed by Drinovec et al. (2015). The MA200 measurements were evaluated based on the instrument's reported status, and flow calibration was performed before each filter change.

Strict lockdown measures were implemented in northern Italy due to the SARS-CoV-2 pandemic from 8 March to 4 May 2020. Therefore, the winter data presented in this study represent a business-as-usual scenario, spanning across two winter seasons. Prior to the analysis, the concentration data were averaged to 1 h to match the time resolution of observed meteorological variables and align with the typical time step used in modelling applications at urban and regional scales.

Absorption measurements from the MA200 devices were further partitioned into specific components, including BC, BrC (hereafter noted as $b_{\mathrm{abs}}^{\mathrm{BC}}$ and $b_{\mathrm{abs}}^{\mathrm{BrC}}$) and their respective sources, namely fossil fuel (FF) and biomass burning combustion (BB), referred to $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$ and $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$, respectively. This partitioning was achieved using the multi-wavelength absorption analyser model (MWAA) developed by Massabò et al. (2015) and Bernardoni et al. (2017b). The MWAA model assumes that the absorption Ångström exponent (α) as defined by Moosmüller et al. (2009), is equivalent for BC and fossil fuel sources while biomass burning is considered the only source of BrC. Under these assumptions, the model holds that the total absorption at each wavelength adheres to the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }{b}_{\mathrm{abs}}=b_{\mathrm{abs}}^{\mathrm{BC}}\left(\mathit{\lambda }\right)+b_{\mathrm{abs}}^{\mathrm{BrC}}\left(\mathit{\lambda }\right)=A\cdot {\mathit{\lambda }}^{-{\mathit{\alpha }}^{\mathrm{BC}}}+B\cdot {\mathit{\lambda }}^{-{\mathit{\alpha }}^{\mathrm{BrC}}}\text{(2)}& {\displaystyle }{b}_{\mathrm{abs}}=b_{\mathrm{abs}}^{\mathrm{FF}}\left(\mathit{\lambda }\right)+b_{\mathrm{abs}}^{\mathrm{BB}}\left(\mathit{\lambda }\right)=A^{\prime }\cdot {\mathit{\lambda }}^{-{\mathit{\alpha }}^{\mathrm{FF}}}+B^{\prime }\cdot {\mathit{\lambda }}^{-{\mathit{\alpha }}^{\mathrm{BB}}}.\end{array}$$

Eqs. (1) and (2) are defined with α^BC=α^FF = 1, which approximates the centre of the α probability density function (PDF) during the morning rush hour at the urban traffic site (Fig. S1 in the Supplement), indicative of fresh, uncoated BC particles, as shown by Liu et al. (2018). The α value for BrC was determined by a preliminary nonlinear fit of Eq. (2), with α^BrC treated as a free parameter (resulting in a mean α^BrC = 3.9). α^BB = 2 was set based on the upper tail of the α PDF calculated from the fit over the five wavelengths and by applying a stringent filter of r² > 0.99 (Fig. S1), as suggested, for example, by Tobler et al. (2021). The same value also aligns with the existing literature data for the Po Valley (Costabile et al., 2017b; Vecchi et al., 2018; Bernardoni et al., 2011). Parameters A and B were then derived for each sample by multi-wavelength fitting of Eq. (1) (with a fixed α^BrC value) and A^′, B^′ by multi-wavelength fitting of Eq. (2). For a comprehensive understanding of the MWAA model, refer to the work of Massabò et al. (2015) and Bernardoni et al. (2017b). In addition, for more detailed information on the source and component assignment performed on the data used in this study, see the companion study by Bigi et al. (2023). Meteorological variables essential for the reconstruction of the wind field by the GRAMM-GRAL model, including downward global radiation and wind speed and direction, were obtained from the meteorological network of the Osservatorio Geofisico di Modena (specifically the CMP, OSS, and POL stations) and from the monitoring network managed by the Regional Agency for Prevention, Environment, and Energy (ARPAE; represented by the DEX station). Table 1 shows the locations and characteristics of the stations used in this study, while Fig. 1 shows their positions on a map. The table also includes information on the locations where soundings were carried out and subsequently used in this study, denoted by the SPC label.

Table 1 Characteristics and locations of meteorological and air quality stations employed in this study.

2.2 The hybrid Eulerian–Lagrangian modelling system

In this study, one of the main objectives was to perform a spatial investigation of BC concentrations in the urban environment of Modena and to provide valuable information on source apportionment. To achieve this, we used a hybrid modelling approach based on two key components: the GRAMM-GRAL dispersion modelling suite, a state-of-the-art Lagrangian simulation system, and the NINFA modelling tool.

The GRAMM-GRAL suite was specifically selected for its ability to simulate the dispersion of emissions within the study area, taking into account the complex urban landscape with its various obstacles such as buildings and structures. On the other hand, NINFA played a crucial role in complementing the output of the Lagrangian model and providing additional insights into BC pollution. NINFA was used to estimate background concentrations from sources outside the urban area and to assess the influence of regional factors on BC levels in Modena. By integrating the strengths of both models, we aimed for a comprehensive understanding of BC concentrations in the study area, capturing the complex dynamics of BC transport in the urban environment. Furthermore, the Lagrangian models used in our approach offered the advantage of generating concentration fields specific to selected source emissions, allowing the estimation of source apportionment and providing valuable information on the contributions of different emission sectors to BC levels. To avoid duplication of emissions in Modena, the BC emission fluxes within the urban domain of interest were set to zero when using NINFA. In addition, BC was treated as an inert species in both models, leading to the exclusion of any chemical reactions involving BC. As a result, the total BC concentration was determined by combining the NINFA output, which represents the contribution of the external areas of Modena, with the GRAL output, which instead represents the contribution of the city. For a more complete understanding of the modelling system and the coupling between an Eulerian and a Lagrangian model used in this study, refer to previously published works such as Veratti et al. (2020, 2021), where a detailed description is provided.

2.3 NINFA set-up

Simulations were conducted using NINFA, a modelling system based on the state-of-the-art chemical transport model CHIMERE v2017r4 (Mailler et al., 2017) and operated by the Regional Agency for Prevention, Environment, and Energy (ARPAE). This suite was employed to investigate surface concentrations of BC aerosol at the regional scale. The study area covered a significant part of the Po Valley, with a grid size of 114 × 86 cells. The grid extended from 43.57–45.79° N to 8.78–13.13° E, with a fine horizontal resolution of 3 km (see Fig. 1 for a visualisation). The vertical grid consisted of 15 levels, ranging from 997 hPa (25 m) for the first layer to 500 hPa. This high-resolution grid was specifically chosen to accurately represent the spatial-scale characteristics of regional air pollution in the area of interest, and it also allowed a precise delineation of the urban area of Modena, where urban emissions were excluded from the simulations. The meteorological input to NINFA was provided by the output of COSMO-2I, a limited-domain atmospheric model used by the Italian National Civil Protection Agency and operated by the local environmental agency ARPAE. The COSMO domain covers the Italian peninsula with a horizontal resolution of 2.2 km and has 30 vertical levels ranging from 20 m to 22 km. Meteorological variables were spatially interpolated to the NINFA grid using the libsim libraries (ARPAE-SIMC, 2023). The aerosol module implemented within CHIMERE provided hourly concentrations of six chemical species including sulfates, nitrates, ammonium, primary and secondary organic particles, dust, and BC. Aerosols were represented using a size bin approach with 10 bins with mean mass median diameters ranging from 10 nm to 40 µm, and BC was simulated with a mass distribution centred at 200 nm with a sigma of 1.2, consistent with previously reported experimental studies (Schwarz et al., 2008; Ning et al., 2013; Li et al., 2023). Detailed information on the gas-phase chemistry and secondary organic aerosol formation scheme used in this study (which is not directly relevant to the ultimate purpose of the current research) can be found in Veratti et al. (2023). In our investigation, we focused only on BC, which was treated as an inert species. Other parameterisations include the van Leer scheme (Van Leer, 1977) for horizontal transport and the Troen and Mahrt formulation (Troen and Mahrt, 1986) for vertical turbulent mixing. In addition, wet and dry deposition processes of aerosols were considered using the equations of Henzing et al. (2006) and Zhang et al. (2001), respectively.

Boundary conditions were obtained from the national model KAIROS (air operational system), which is the air quality forecast model operated by the regional air quality agency ARPAE on a daily basis throughout Italy (Stortini et al., 2020), using the same version of CHIMERE used in this study with a horizontal spatial resolution of 7 km. The simulation period chosen for our study was from 4 February to 7 March 2020 and from 26 December 2020 to 21 January 2021, corresponding to the availability of BC measurements in Modena. In order to ensure the reliability and accuracy of our simulations, we implemented a spin-up period of 3 d in order to reach an equilibrium state and to establish stable initial conditions before the periods of interest.

In order to perform a source-tagged BC simulation aimed at distinguishing the contribution to total BC concentrations from different sources, we made modifications to the original CHIMERE code. These modifications allowed the specific allocation of the contributions of different sectors coming from outside the urban area of Modena.

The performance of CHIMERE in reproducing gaseous species, total particulate matter (PM), and PM components has been extensively evaluated in several intercomparison studies. Notably, the EURODELTAIII (Bessagnet et al., 2016; Mircea et al., 2019), AQMEII (Pirovano et al., 2012), and POMI (Pernigotti et al., 2013) exercises have been significant in this context. Specifically, Mircea et al. (2019) focused on the simulation of carbonaceous aerosols over Europe. When EC concentrations within the PM_2.5 matrix were compared with corresponding measurements, CHIMERE showed satisfactory performance in reproducing the observed trend for the years 2006–2009, with a Pearson correlation coefficient (r) between modelled and measured concentrations ranging from 0.6 to 0.85, a normalised centred root-mean-square error (RMSE) between 1.5 and 2.0, and a normalised standard deviation (SD) between 0.05 and 1.25.

In the EURODELTAIII exercise (2006–2009), CHIMERE generally overestimated O₃ concentrations compared to other participating models, with a mean bias (MB) between 6.3 and 22.5 µg m⁻³. However, in 2009, all models, including CHIMERE, underestimated the measured concentrations due to biased boundary conditions. Despite this, CHIMERE had the lowest correlation coefficient (ranging from 0.27 to 0.71) and also the lowest RMSE among the models. For NO₂, its performance was comparable to other models, with r ranging from 0.67 to 0.72 and MB ranging from −0.64 to 0.64 µg m⁻³, depending on the simulation year. For SO₂, the correlation coefficient ranged from 0.2 to 0.4, similar to other models, but CHIMERE was closest to the observations (MB from −0.46 to 0.13 µg m⁻³). For PM₁₀, CHIMERE generally underestimated the measurements, as did other models, with MB ranging from −6.59 to −0.52 µg m⁻³, but it achieved the highest correlation coefficient (0.7). Similarly, CHIMERE performed well for PM_2.5 over all simulated years, with the lowest RMSE (6.59 µg m⁻³) and a comparable MB (between −1.00 and −2.39 µg m⁻³) to the other models.

In addition, in the POMI intercomparison exercise, CHIMERE was used alongside EMEP, AURORA, CAMx, TCAM, and REM-CALGRID to simulate O₃ and PM over the Po Valley. The results indicated that CHIMERE's performance was comparable to the other models. For daily PM₁₀, CHIMERE had one of the lowest RMSE values (29.1 µg m⁻³) and the highest correlation coefficient (0.6), with other statistical metrics in good agreement with the other models. However, for O₃, CHIMERE had the highest absolute MB (−9.1 µg m⁻³) and also one of the lowest RMSE values (31.5 µg m⁻³), with other metrics in agreement with the other models.

2.4 GRAMM-GRAL description and set-up

The GRAMM-GRAL modelling system is a comprehensive approach designed for simulating air pollution in urban areas. It integrates two tools: GRAMM (the Graz Mesoscale Model) and GRAL (the Graz Lagrangian Model). While GRAMM focuses on the simulation of meteorological variables such as temperature, humidity, wind speed, and wind direction, GRAL is a Lagrangian particle dispersion model capable of reproducing pollutant concentrations in the atmosphere. This system is particularly well suited for simulating air pollution in urban environments, even in the presence of steep topography.

GRAMM is a non-hydrostatic model that addresses the conservation equations for mass, enthalpy, momentum, and humidity. It takes into account the effects of topography and surface interactions, such as the exchange of heat, momentum, humidity, and radiation, across different land use categories. To handle turbulence, an algebraic turbulence model (Pandolfo, 1969) is used when gradient Richardson numbers exceed 0.3, while a k–ϵ approach is used for lower values. Standard vertical profiles of wind and temperature are used to initiate and drive GRAMM at its boundaries, depending on the synoptic-scale forcing and stability class chosen as input. The initialisation process follows the Pasquill–Gifford classification, and a detailed theoretical formulation and a practical implementation can be found in the works of Almbauer et al. (1995) and Oettl (2015b). The ability of GRAMM to reproduce meteorological conditions with different topographic settings and resolutions has been demonstrated in previous studies such as Almbauer et al. (1995), Oettl (2015b), Thunis et al. (2003), Oettl (2021), and Oettl and Reifeltshammer (2023). Furthermore, its performance has recently been evaluated by comparison with the widely used Weather Research and Forecasting (WRF) model (Skamarock et al., 2008), showing very similar results for wind and temperature fields compared to observations (Oettl and Veratti, 2021).

Within GRAMM, the microscale model GRAL is embedded through a one-way coupling from the larger to the smaller scale. GRAL offers the option to operate in either diagnostic or prognostic mode. In the diagnostic mode, the flow field around buildings is determined by interpolating the large-scale wind fields from GRAMM output assuming a logarithmic wind profile close to urban surfaces; then, a Poisson equation is applied to correct velocities and guarantee the mass conservation. On the other hand, the prognostic mode, despite its higher computational cost, ensures superior accuracy by explicitly computing the flow through forward integration of the Reynolds-averaged Navier–Stokes equations (RANS; Oettl, 2015c), while turbulence is estimated using a standard k–ϵ approach, neglecting Coriolis and buoyancy forces (Oettl, 2015c).

In this research, a cascade of scales was used to determine wind patterns within the city, from the synoptic to the building level. The study started with the simulation of large-scale meteorological conditions in the region surrounding the city of interest. An area of 30 km × 30 km was considered, centred on Modena, with a resolution of 200 m (see Fig. 1). The simulation was carried out using GRAMM version 22.09, taking into account the local topography obtained from the Geoportale-Emilia-Romagna (2023) and land use data extracted from the Corine land cover database updated to 2018 (CCL, 2018). These mesoscale flow patterns were used as the driving force for the computation of high-resolution winds within the city, specifically accounting for buildings and street canyons. Finally, Lagrangian dispersion calculations were performed using the dispersion module of GRAL (version 22.09), with constraints imposed by the previously generated microscale wind fields. The GRAL domain includes the historical city centre, the ring road, and most of the urban environment of Modena, covering a domain of approximately 7.9 km × 6.5 km, with a horizontal resolution of 4 m (see Fig. 1).

Since computation of the city-scale wind field at a 4 m resolution over an entire urban environment is computationally intensive, the approach of Berchet et al. (2017) was implemented using a discrete representation of different weather situations. This approach took into account the classification of atmospheric stability, large-scale wind speed, and direction at the boundaries of the domain into discrete categories. Large-scale wind directions were divided into 36 classes of 10° each, while wind speeds were classified into 15 classes ranging from 0.25 to 7 m s⁻¹. Atmospheric stability was classified according to the Pasquill–Gifford classes (EPA, 2000). This resulted in a catalogue of about 1100 physically meaningful reference weather situations that occurred during the two periods of interest. Once the catalogue of meteorological situations was completed, at each simulation time step, the wind field was selected by matching the meteorological catalogue with the observed meteorological conditions at the reference urban stations within the urban area of Modena (see Fig. 1). Lagrangian dispersion calculations were then performed at hourly intervals in prognostic mode, using a particle diameter of 2.5 µm. In order to validate the model wind fields at the urban scale, we used the same stations that were used to select the wind field situation from the catalogue. The location of the stations is shown in Fig. 1, while their characteristics are reported in Table 1. On the other hand, the validation of the results obtained with NINFA and GRAL was carried out using eBC measurements at two sites, the urban background and traffic sites, respectively.

2.5 Anthropogenic emissions

2.5.1 City scale

The anthropogenic emissions used in this study were obtained from the regional emission inventory for the Emilia–Romagna region, regularly compiled by the local environmental agency ARPAE (INEMAR, 2023). This inventory provides annual totals of PM₁₀, PM_2.5, NO_x, CO, SO₂, NH₃, and non-methane volatile organic compounds (NMVOCs) for each municipality, in this case specifically for the year 2017. For transport emissions, however, we used a bottom-up approach rather than relying solely on the regional inventory. The latter has been shown to be successful in reporting high-resolution traffic emissions in previous studies conducted in the same region (Ghermandi et al., 2019, 2020; Veratti et al., 2020, 2021).

To estimate traffic emissions at the city scale, we employed a comprehensive bottom-up approach that integrates emission factors (EFs) and activity data. The activity data were derived from traffic flows simulated by the PTV VISUM model. Then, specific exhaust EFs were applied to estimate the total emissions, taking into account various factors such as fleet composition, vehicle type, fuel, engine capacity, load displacement, road slope, Euro emission standard, and average travelling speed. Additionally, we considered non-exhaust traffic EFs, which are contingent upon vehicle weight and travelling speed, to estimate emissions from tyre, brake, and road wear.

To estimate BC exhaust emissions from road transport, we adopted the tier 3 methodology outlined in the European Monitoring and Evaluation Programme/European Environment Agency (EMEP/EEA) guidelines (Ntziachristos and Samaras, 2019), which includes the calculation of both hot emissions, which occur when the engine is operating at its normal temperature, and emissions during transient thermal engine operation, commonly referred to as cold start emissions. To determine the EF for each vehicle fleet category, we used flow velocity estimates from PTV VISUM, and the individual EFs were averaged taking into account the local fleet composition (ACI, 2023) and the total annual kilometres travelled by each vehicle category considered (ISPRA, 2023). As a result, we obtained two EFs per road segment: one for light vehicles (such as cars and mopeds) and another for light and heavy commercial vehicles. This estimate corresponds to the two classes available in the PTV VISUM output and allows a consistent representation of BC exhaust emissions.

On the other hand, for the estimation of non-exhaust BC traffic emissions, the tier 2 approach described in Ntziachristos and Boulter (2019) was implemented. The calculation was based on the so-called detailed methodology, which takes into account the speed dependency of tyre and brake wear emissions from moving vehicles, while for road surface wear, the emissions depend only on the number of vehicles travelling on a given road. The exhaust and non-exhaust calculation methods used in this study were coded in the R programming language and embedded in a package called VERT (Vehicular Emissions from Road Traffic; Veratti et al., 2024). The VERT package provides a comprehensive framework for estimating emissions from road traffic, including exhaust, non-exhaust and evaporative emissions. The calculation process takes into account fleet composition and vehicle flows for cars, mopeds, light-duty vehicles, and trucks or aggregated classes depending on the level of data available. This package includes the latest emission factors and methodologies recommended by the European EMEP/EEA guidelines for not only BC emissions but also a range of other pollutants including CO, NO, primary NO₂, NO_x, VOC, CH₄, PM, Organic Carbon, N₂O, NH₃, CO₂, and SO₂. Furthermore, in recognition of the significant uncertainties associated with PM speciation in the determination of BC emission factors, VERT provides an uncertainty range for both exhaust and non-exhaust emission determination based on the speciation interval specified in the methodology documentation (Ntziachristos and Samaras, 2019; Ntziachristos and Boulter, 2019). In this study, we aimed to assess the impact of emission uncertainties on the model results and to investigate the potential range of outcomes. To achieve this, we ran a series of scenarios incorporating different hypotheses regarding BC speciation factors and the distribution of emission sources within the city boundaries. This approach allowed us to obtain a comprehensive set of BC concentration maps covering a plausible range of outcomes.

For traffic emissions, we ran a base case scenario using the reference speciation factor for BC as suggested by the European EMEP/EEA guidelines and two additional scenarios using the lower- and upper-range values for BC speciation. In addition, to increase the comprehensiveness of our assessment of the overall impact of traffic, we extended the emission estimates beyond exhaust and non-exhaust BC emissions and added resuspension due to traffic circulation. To include this source, we conducted five additional scenarios. Four of these were derived from the methodology proposed by Amato et al. (2012) and implemented as described in the HERMESv3 model (Guevara et al., 2020), while the fifth was evaluated based on the formulation of the NORTRIP model by Denby et al. (2013). The Amato et al. (2012) approach relies on vehicle-specific emission factors, the values of which are detailed in Table 2. On the other hand, the NORTRIP formulation is based on ensuring mass conservation on the road surface. For this reason, we used a three-step procedure to estimate the resuspension effects using the NORTRIP model. First, we calculated the dispersion of the full set of emissions, including traffic exhaust, non-exhaust, and emissions from other activity sectors. Then we used the BC deposition on the road surface obtained from this calculation to estimate the resuspension effects based on the formulation reported below (Eq. 3). Finally, we performed another simulation using these emission rates to determine the contribution of resuspension to atmospheric concentrations. The calculation of the resuspension rate Q_resusp, which represents the number of particles resuspended into the atmosphere according to the NORTRIP model, is determined according to Eq. (4), taking into account the mass deposited on the road surface M_dep.

$$\begin{array}{}\text{(3)}& Q_{\mathrm{resusp}}=M_{\mathrm{dep}}\cdot F_{\mathrm{resusp}},\end{array}$$

where F_resusp is the resuspension factor and can be calculated as

$$\begin{array}{}\text{(4)}& F_{\mathrm{resusp}}=\sum _{v=\mathrm{1}}^{\mathrm{2}}{N}_{\mathrm{v}}\left({\displaystyle \frac{u_{\mathrm{v}}}{u_{\mathrm{ref}\left(r\right)}}}\right)f_{\mathrm{0},\mathrm{v}}.\end{array}$$

The variables in Eq. (4) are defined as follows: v represents the vehicle type, N_v represents the vehicle flow (measured in vehicles per hour), u_v represents the vehicle speed (measured in km h⁻¹), u_ref(r) represents the reference vehicle speed for the resuspension process (measured in km h⁻¹), and f_0,v represents the reference mass fraction of the resuspension process per vehicle. In line with previous studies by Thouron et al. (2018), Denby et al. (2013), and Lugon et al. (2021), this study adopts the values u_ref(r) = 50 km h⁻¹, f_0,HDV = $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ per vehicle, and f_0,LDV = $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ per vehicle. All traffic scenarios simulated in this study are summarised in Table 2.

Table 2 Simulated emission scenarios for traffic sources.

^a See Table 3-91 of Ntziachristos and Samaras (2019) for the exhaust PM speciation factors for different vehicle technologies.
^b See Tables 3-4 and 3-6 of Ntziachristos and Boulter (2019) for the non-exhaust emission factor range for different vehicle categories.

A combination of speciation factors and detailed inventory data was used to estimate emissions from other activity sectors. As introduced at the beginning of this section, the ARPAE inventory is used as a reference, while activity-dependent BC speciation factors reported by the EEA (2019) and listed in Table S2 were used to convert total PM emissions into BC estimates. In order to distribute the emissions within the urban area of Modena, land use data and building characteristics were used as proxy variables. However, it should be noted that in the Emilia–Romagna region there are still significant uncertainties regarding the distribution of emissions from biomass burning for non-industrial combustion. To address this limitation, we developed five different emission scenarios to capture the variability in emissions in terms of spatial distribution. These scenarios considered different source locations and different total emissions allocated to different areas of the city, including the historic centre, discontinuous urban fabric, and rural areas. Table 3 summarises the characteristics of each scenario.

Table 3 Simulated emission scenarios for non-industrial combustion sources.

The daily temporal modulation of non-industrial combustion emissions is based on the concept of heating degree days, a metric developed to capture the energy demand required to heat a structure based on both external and internal temperatures. The formula used in this study to incorporate this calculation is derived from Mues et al. (2014). In addition to transport and biomass combustion, other relevant SNAP sectors included in the calculation are emissions from industry and other mobile machinery. However, due to their relatively small contribution in annual totals compared to transport and biomass burning, no additional scenarios were simulated to assess their impact range.

2.5.2 Regional scale

The emission data used for NINFA were taken from the same inventory used at the urban scale, where the total annual emissions for each pollutant are reported at the municipal level. To ensure a precise depiction at a more detailed spatial and temporal resolution, a downscaling procedure was implemented. Emissions were spatially allocated at the model resolution using the Corine land cover database (CCL, 2018) and temporally distributed using monthly, daily, and hourly profiles typical of northern Italy (Veratti et al., 2023). To avoid double counting BC emissions within the urban area of Modena, the emission fluxes from this area were set to zero.

The availability of high-resolution emission inventories produced by environmental agencies can sometimes cause challenges. Therefore, a comparison between the annual BC emissions used in a study for the Emilia–Romagna region (INEMAR 2017) and common European and global emission datasets can provide valuable insights for future modelling studies in the Po Valley. Table 4 presents a comparison between INEMAR 2017 (INEMAR, 2023; Marongiu et al., 2024), CAMS-REGv4.2 (2018; Kuenen et al., 2022), EDGARv6.1 (2018; Crippa et al., 2020), and EMEP (2017; Ullrich et al., 2023), categorised by sector.

The comparison results show that for non-industrial combustion, EDGARv6.1 and EMEP agree quite well with INEMAR, with discrepancies between 8 % and 10 %, while CAMS-REGv4.2 is about 60 % higher than INEMAR. For traffic emissions, the annual totals are very close to each other, with differences ranging between 3 % and 22 %, confirming that this sector is relatively well constrained. However, total emissions for other sectors exhibit significant discrepancies, partly due to the different classification systems used by different inventories to categorise emissions. Despite these differences, total BC emissions are relatively well aligned, with differences ranging from 11 % to 17 %.

Table 4 Comparison of annual BC emissions for the Emilia–Romagna region across different emission sectors as reported by INEMAR 2017, CAMS-REGv4.2, EDGARv6.1, and EMEP 2017.

To assess the influence of external areas on BC concentrations in Modena, six scenarios were developed. Three scenarios focused on the speciation factors of both exhaust and non-exhaust traffic emissions. In particular, to convert the total PM emissions reported in the inventory into BC emissions, vehicle-fleet-dependent BC speciation factors calculated with the VERT package were used. The first scenario used average BC conversion factors, while the other two scenarios represented the lower and upper bounds of the BC speciation values. For non-industrial-combustion BC emissions, three different spatial allocations were considered. In the first allocation, municipal totals were distributed with 40 % allocated to rural areas and 60 % to urban areas. In the second allocation, the same proportion was maintained, but 60 % was allocated to rural areas and 40 % to urban areas. In the third allocation, 80 % was allocated to rural areas and 20 % to urban areas. Table 3 provides an overview of each scenario simulated at the regional level.

3 Model evaluation

A number of statistical metrics have been used to evaluate the hybrid modelling system against observed data. These metrics include MB; normalised mean bias (NMB); r; the proportion of predicted values within a factor of 2 of the observations, called factor of 2 (FAC2); the normalised mean square error (NMSE); fractional bias (FB); the normalised average difference (NAD); and RMSE. Detailed definitions of these metrics can be found in Sect. S3 of the Supplement.

4 Results

4.1 Measurements from MA200

Figure 2 shows the average diurnal variability in the aerosol absorption at 880 nm (panel a) from the two MA200 instruments partitioned using the results of the MWAA model, distinguishing between fossil fuels, biomass combustion, and BrC. The same partitioning was also used to differentiate eBC concentrations (Fig. 2b) between fossil fuels (eBC-FF) and biomass combustion (eBC-BB), applying a fixed MAE (see Sect. 2.1). Due to the lack of in situ validation measurements to assess the MAE of BrC and the significant uncertainty associated with its determination (see, for example, the diverse structures and properties of the wide range of chemical organic compounds coupled with its sensitivity to variations in time and location), we chose not to convert $b_{\mathrm{abs}}^{\mathrm{BrC}}$ into concentrations. This decision was taken to avoid introducing additional uncertainty associated with assuming a specific MAE for BrC.

Figure 2 The average daily cycle of $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$, and $b_{\mathrm{abs}}^{\mathrm{BrC}}$ is shown for panel (a), while panel (b) displays eBC-BB and eBC–FF concentrations, all apportioned using the MWAA model. The solid lines represent the average absorptions or concentrations, while the shaded areas indicate the interquartile range. The top graphs illustrate concentrations measured at the urban background site, while the bottom graphs display concentrations measured at the traffic site.

The measurements reported in Fig. 2 were carried out during two different periods. The first period lasted from 4 February to 7 March 2020, while the second period spanned from 26 December 2020 to 21 January 2021 for the traffic site and from 26 December 2020 to 7 January 2021 for the urban background site.

When analysing the results, it can be seen that $b_{\mathrm{abs}}^{\mathrm{BrC}}$ remain consistently very low on average throughout the day for both sites, averaging around 0.5 Mm⁻¹ for the first period and around 1.2 Mm⁻¹ for the second period. On the other hand, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$ shows the characteristic behaviour associated with traffic emissions, with two peaks during rush hours. Specifically, these peaks occur in the morning between 07:00 and 08:00 GMT and in the evening between 18:00 and 19:00 GMT or slightly delayed at the urban background site for the first period, between 18:00 and 20:00 GMT. This pattern is also observed for other traffic-related tracers such as NO, NO₂, and benzene measured at the same two sites, not shown here but reported in a companion study by Bigi et al. (2023).

In contrast, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ shows less variability and generally lower diurnal absorptions compared to $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$. In particular, their behaviour is different, with $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ reaching a minimum in the early afternoon, typically between 15:00 and 17:00 GMT. This phenomenon can be attributed to a combination of factors, including a higher mixing layer depth and lower emissions during these hours. In the later part of the afternoon, after 17:00 GMT, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ typically increase and reach a maximum between 20:00 and 21:00 GMT, when the emitting sources are likely to be at their peak emission levels and the mixing layer depth is gradually decreasing due to surface cooling. Despite a likely decrease in the intensity of biomass burning emissions during the night, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ remain relatively high due to a temperature inversion that often characterises the area, while during the rest of the day $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ gradually decrease until reaching the minimum explained above.

Comparing the first and second periods, it is interesting to note that the mean absorption coefficient for BC in the second period was significantly higher than that of the first, with mean values for $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$ (eBC-FF) and $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ (eBC-BB) of 19.7 ± 15.1 Mm⁻¹ (1.94 ± 1.49 µg m⁻³) and 4.5 ± 5.9 Mm⁻¹ (0.44 ± 0.58 µg m⁻³), respectively, at the traffic site for the first period and 29.4 ± 24.6 Mm⁻¹ (2.90 ± 2.43 µg m⁻³) and 13.3 ± 11.2 Mm⁻¹ (1.32 ± 1.11 µg m⁻³) for the second period at the same station. $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$ (eBC-BB) at the urban background site expresses similar behaviour with a mean value equal to 5.2 ± 4.7 Mm⁻¹ (0.52 ± 0.46 µg m⁻³) for the first period and to 9.6 ± 8.7 Mm⁻¹ (0.95 ± 0.86 µg m⁻³) for the second period. A possible explanation for this increase can be found in the different meteorological conditions in the two periods. More specifically, the second period shows more-stagnant conditions compared to the first period. This is due to more stable atmospheric conditions, characterised by recurrent thermal inversions (Figs. S9 and S10) and a higher frequency of stable atmospheric classes compared to the first period (Figs. S2 and S3), which likely led to lower-mixing-layer height. The frequency of occurrence of stable atmospheric classes (G and F, according to the Pasquill and Gifford formulation) is 0.53 for the first period and 0.63 for the second period. In addition, the mean hourly temperature observed at the OSS station, representing the urban area, was significantly lower in the second period (2.9 °C) compared to the first period (9.4 °C). This probably led to an increase in biomass burning for domestic heating, contributing significantly to the total eBC concentrations (see Fig. 2).

On the other hand, $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$ (eBC-FF) at the urban background site shows a different trend. Contrary to the other site and $b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$, its mean value decreased from 9.8 ± 7.1 Mm⁻¹ (0.97 ± 0.70 µg m⁻³) in the first period to 7.3 ± 5.5 Mm⁻¹ (0.72 ± 0.55 µg m⁻³) in the second period. The most plausible explanation for this discrepancy lies in the type of days included in the second period. Specifically, the measurements at the urban background site were taken over a period of 13 d, 7 of which coincided with holidays during Christmas time, resulting in lower traffic volumes throughout the city. It is therefore reasonable to expect a limited contribution from traffic during this period. Furthermore, it is worth noting the impact of strong stagnation conditions that occurred in Modena between 13–14 and 18–19 January 2021 (see Figs. S2 and 5), verified after the end of the measurements at the urban background site.

Numerous studies have reported eBC levels across Europe. In recent years, for example, Savadkoohi et al. (2023) carried out a comprehensive comparison of eBC levels at 50 monitoring sites, including the cities of Paris, Milan, Bern, and Lille for the periods of 2016–2019, 2019–2021, 2015–2021, and 2017–2019, respectively. At traffic sites in Paris and Milan, the eBC-FF concentrations (95 % confidence interval) were 2.17 (0.55, 5.06) and 2.68 (0.75, 6.11) µg m⁻³, while the eBC-BB concentrations at the same two sites were 0.31 (0.00, 0.99) and 0.53 (0.07, 1.88) µg m⁻³, respectively. At urban background sites, Bern and Lille had eBC-FF concentrations of 0.78 (0.17, 1.93) µg m⁻³ and 0.68 (0.14, 1.76) µg m⁻³, while eBC-BB concentrations were 0.28 (0.04, 0.79) and 0.33 (0.04, 1.01) µg m⁻³, respectively. In addition, Pashneva et al. (2024) reported eBC-FF and eBC-BB concentrations of 0.78 ± 0.76 and 0.35 ± 0.42 µg m⁻³ for the city of Vilnius during winter 2021–2022. On the other hand, other urban locations in Europe such as Athens, Madrid, and Rome showed average total eBC values of 2.8, 2.3, and 2.9 µg m⁻³, respectively, for the winter periods between 2015 and 2019, between 2014 and 2015, and for January and February 2017 (Liakakou et al., 2020; Becerril-Valle et al., 2017; Costabile et al., 2017a).

Table 5 shows the average results of the MWAA model, providing insight into the partitioning of components (BC vs. BrC) and sources (FF vs. BB) for both the first and second periods of our investigation. Furthermore, the table shows the results of the MWAA and aethalometer model partitioning, as well as the results of the application of the EPA PMF 5.0, using hourly absorption data and elemental concentrations, as observed in previous research carried out in the Po Valley. The comparative analysis underlines the consistency of the results obtained in our study using MA200 instruments, which are on average in line with the findings of other studies carried out in the same region.

Furthermore, source apportionment results in Modena are consistent with those from other European urban areas. For instance, Savadkoohi et al. (2023) reported FF and BB partitioning comparable to Modena at traffic sites in Stockholm, Paris, and Milan during winter. For the periods of 2014–2019, 2016–2019, and 2019–2021, the FF and BB partitioning performed with the aethalometer model were 80 %$/$20 %, 81 %$/$19 %, and 69 %$/$31 %, respectively. Similarly, a traffic site in Bern showed comparable source apportionment to Modena during the winter season in 2014–2018, with an FF and BB partitioning of 60 %$/$39 % (Grange et al., 2020). For urban background locations, Savadkoohi et al. (2023) also presented results for the cities of Zurich, Athens, Paris, and Bucharest, which were similar to those observed for the urban background location in Modena. The FF and BB partitioning during winter for these cities was 67 %$/$31 %, 59 %$/$41 %, 56 %$/$44 %, and 50 %$/$49 % for the periods of 2012–2021, 2017–2019, 2016–2019, and 2014–2022, respectively. Other studies reported similar percentages for urban background locations, such as Vilnius in winter 2021–2022 (68 %$/$31 %; Pashneva et al., 2024), the surroundings of Šibenik (Croatia) in February 2019 (68 %$/$31 %; Milinković et al., 2021), the urban background of Ruhr (Germany) in winter 2013–2014 (50 %$/$50 %), and Helsinki in winter 2015–2017 (54 %$/$46 %; Helin et al., 2018).

Bernardoni et al. (2021)Mousavi et al. (2019)Forello et al. (2019)Ferrero et al. (2018)

Table 5 Summary table of component (BC vs. BrC) and source (FF vs. BB) apportionment based on this work and other studies in the literature conducted in the Po Valley. Most values are reported at 880 nm, with the exception of Forello et al. (2019), who present data at 780 nm. RB, UB, and UT stand for rural background, urban background, and urban traffic, respectively.

^a Seven wavelengths (370, 470, 520, 590, 660, 880, and 950 nm) fit with fixed multiple-scattering enhancement parameters.
^b Four wavelengths (470, 520, 660, and 880 nm) fit with fixed multiple-scattering enhancement parameters.
^c Five wavelengths (470, 520, 590, 660, and 880 nm) fit with fixed multiple-scattering enhancement parameters.
^d Source apportionment was conducted using EPA PMF 5.0, incorporating hourly absorption data at four distinct wavelengths (405, 532, 635, and 780 nm), in combination with hourly elemental concentrations from samples gathered and analysed using the same filters. The outcomes of this analysis revealed that the third factor in the apportionment could be attributed to sulfate, constituting approximately 20 % of the overall composition. The fourth factor, on the other hand, was primarily composed of resuspended dust, accounting for approximately 4 % of the total, with other factors of lesser significance contributing approximately 1 %.
^e Aethalometer model.
^f MWAA model.
^g Total BC, representing the cumulative sum of BC-FF and BC-BB.

4.2 Modelling results

4.2.1 Meteorology

Using the GRAMM-GRAL modelling system, a comprehensive catalogue of mesoscale and microscale flow patterns was calculated. Following the methodology outlined in Sect. 2.4, a series of hourly wind fields were generated to drive the Lagrangian dispersion in Modena. Whenever possible, meteorological observations from the stations shown in Fig. 1 were used. First, the wind fields were computed independently of any observation, following the discretisation procedure described in Sect. 2.4. Subsequently, the wind fields were selected to closely match the actual conditions. It is important to acknowledge that although a given wind field from the catalogue may not be perfectly consistent with all stations simultaneously, this study relies exclusively on data from stations located in the urban area of Modena to construct the catalogue. This is due to the unavailability of measurements from other locations within the GRAMM domain and within a range of 50–60 km. Furthermore, considering that the area of interest is almost flat, we expect that this set-up will not significantly affect the representation of the flow fields in the area of interest. Figure 3 shows the comparison between modelled and measured hourly wind speed for the available urban stations for the period of interest. Panel (a) represents the first period, while panel (b) shows the wind speed comparison for the second period. Despite some notable underestimations of the observed wind speed on 5, 10, 26, 27, and 28 February 2020 due to the wind category discretisation process, GRAMM-GRAL generally reproduced the trend quite well for both periods. For the first period, the wind speed is generally underestimated, with MB between −0.44 and −0.27 m s⁻¹ at CMP and OSS, respectively (corresponding to −34 % and −9 % of the NMB), while at DEX, the NMB is approximately zero and the MB is 0.02  m s⁻¹. Similarly, during the second period, the wind speed is generally underestimated at CMP and POL, with MB values of −0.28 and −0.03 m s⁻¹, respectively, while it is overestimated at OSS with an MB of 0.31 m s⁻¹. The corresponding NMB values are −37 %, −2 %, and +20 %, respectively. The RMSEs of the wind speed simulated by GRAMM-GRAL are 0.88, 0.72, and 1.16 m s⁻¹ for CMP, OSS, and DEX during the first period and 0.60, 0.83, and 0.89 m s⁻¹ for CMP, POL, and OSS during the second period. These RMSE values are in line with the recommended statistical benchmark suggested by the EEA guidelines for meteorological wind field assessment (EEA, 2011), which suggests RMSE values of less than 2 m s⁻¹. Other statistical indices, including FAC2, NMSE, and r, are reported in Table 6, and the related scores are consistent with similar studies conducted in the same area during previous simulation years (Veratti et al., 2020, 2021). To complement the wind speed analysis, Fig. S4 displays wind roses comparing the modelled and observed wind speed and direction for the same stations discussed above. The wind roses confirm that the modelling system effectively reproduces many of the observed features, with the observed winds being generally captured well by GRAMM-GRAL, especially during the second period.

Figure 3 Hourly time series of measured wind speed (black) at four urban meteorological stations in Modena alongside microscale modelled values by GRAMM-GRAL (red) for two periods: (a) 4 February to 7 March 2020 and (b) 26 December 2020 to 21 January 2021.

Table 6 Statistical analysis of hourly wind speed computed at urban meteorological stations using the GRAMM-GRAL modelling system.

Further analysis was carried out on the height of the planetary boundary layer (PBL), a key driver of atmospheric dispersion. Although direct measurements of the PBL height were not available for Modena, we made both quantitative and qualitative comparisons with observations taken in rural areas at the S. Pietro Capofiume (SPC) station, approximately 50 km east of Modena (see Fig. 1 and Table 1). Specifically, we estimated the PBL height using the Richardson number (Ri) derived from sounding data at 00:00 and 12:00 GMT, using 0.25 as the critical value for identifying turbulent conditions (Lyons et al., 1964; Galperin et al., 2007; Grachev et al., 2013). These estimates were then compared with the PBL height simulated by NINFA, which also uses the Richardson number in its calculations (Troen and Mahrt, 1986). Figures S5 and S6 provide an overview of this comparison for two periods: 15 February to 7 March 2020 and 26 December 2020 to 21 January 2021. For the first period, data were only available at 00:00 GMT, while for the second, data were also available at 12:00 GMT. The results show that between 15 February and 7 March 2020, NINFA generally underestimates the sounding estimates, with an MB of −100 m (−52 %), resulting in a shallower mixing layer compared to the measurements. However, during the second period, the PBL height was reproduced better, with a limited MB of −47 m (−27 %) at 00:00 GMT and −38 m (−12 %) at 12:00 GMT, showing robust performance in simulating the vertical structure of the atmosphere. Despite the more pronounced underestimation of the first period, NINFA showed similar performance to other meteorological models applied in the Po Valley and other locations in Italy in reproducing the PBL height derived from soundings (Ferrero et al., 2011; Avolio et al., 2017).

For a qualitative comparison, the PBL height simulated by GRAL over the urban area of Modena is also included in Figs. S5 and S6. Although a quantitative analysis between GRAL and soundings is not possible due to the fact that the measured area is outside the GRAL domain and various factors may cause differences in the vertical turbulence profile between urban and rural areas (e.g. urban heat island effects and anthropogenic heat sources), this comparison serves as a basis for hypothesising whether GRAL can realistically reproduce the PBL height during sounding time.

The results of the GRAL simulations show that during the first period, the PBL height is in agreement with the sounding data on most days (16 out of 22). However, on 6 of the 22 d, GRAL significantly overestimates the measurements at 00:00 GMT, with values up to 800 m, which corresponds to the domain top internally set in the model code. During the second period, the PBL height modelled by GRAL generally matches the sounding data at both 00:00 GMT and 12:00 GMT, except for 9 d when there is a significant overestimation (up to 4 times the observations) at 12:00 GMT. A more detailed analysis of these episodes shows that during the night (00:00 GMT) when stable or very stable atmospheric conditions are imposed on GRAL, the PBL height tends to match the values observed at SPC. Conversely, under neutral conditions, the PBL height simulated by GRAL overestimates the observations by a factor of 3 to 6. This overestimation under neutral conditions is also evident at 12:00 GMT during the second period, with 6 out of 9 episodes occurring under these conditions.

Despite the challenges faced by GRAL, a common limitation for urban air quality models (Sokhi et al., 2022), the simulated PBL height appears realistic on most days. In addition, it is important to note that the limited number of observational data points prevents further conclusions from being drawn for other times of the day.

4.2.2 Comparative analysis of modelled BC concentrations using a hybrid Eulerian–Lagrangian modelling system in Modena

The background BC concentrations derived from NINFA in the urban area of Modena were integrated with the BC concentrations simulated by the GRAMM-GRAL modelling system, which specifically considers emissions only at the urban scale. The scenarios considered for the evaluation include TRF-exh-1 for traffic exhaust emissions and TRF-nexh-1 for non-exhaust emissions. In the case of resuspension and biomass burning, which are characterised by a high level of uncertainty, the average concentrations calculated from all simulated scenarios were taken into account. For NINFA, the scenarios TRF-bck-1 and DMH-bck-1 were considered. This specific emission configuration will be referred to as the base case scenario.

To evaluate the effectiveness of the hybrid modelling system in various urban settings, modelled concentrations at station locations were compared with observed eBC levels from MA200 devices. Figures 4 and 5 show the time series comparing measured (eBC) and modelled (BC) concentrations, divided into their contributors (biomass burning, fossil fuels, and the sum of the two) for the first and second periods. The analysis reveals notably robust performance of the hybrid system in modelling concentrations during the first simulated period (5 February–8 March, 2020), despite BB generally being overestimated at the traffic site. The bias of total modelled BC concentrations is very low at both sites: MB is −0.12 µg m⁻³ for urban traffic (corresponding to −5 % of NMB) and 0.03 µg m⁻³ for urban background (corresponding to +2 %). The Pearson correlation coefficient between modelled and observed concentrations is 0.62 at the traffic site and 0.51 at the urban background station. In addition, the modelling system generally captured the daily peaks correctly, even under different meteorological conditions. In particular, several relatively high wind speed episodes for typical Modena meteorological conditions occurred during this period, such as on 5, 11, 12, 26, and 28 February and on 2 and 5 March 2020, with speeds reaching 11.2 m s⁻¹. The modelling system showed adaptability in simulating both calm and windy periods. Conversely, performance of the model in reproducing observations decreased during the second period (26 December 2020–21 January 2021), as observed eBC concentrations experienced significantly higher values than in the first period. The Pearson correlation coefficients between modelled and observed concentrations are 0.34 and 0.38 at the traffic and background sites, respectively. The corresponding MBs are −0.92 µg m⁻³ (corresponding to −22 % of NMB) and +0.25 µg m⁻³ (corresponding to +15 % of NMB) at the same two stations. Favourable meteorological conditions for the accumulation of pollutants were observed in both Modena and the surrounding areas. These conditions were mainly characterised by high-pressure systems and persistent thermal inversions in the lower atmospheric layers, resulting in a likely shallow mixing layer height.

Figure 4 Hourly time series of observed (eBC) and simulated (BC) concentrations at urban traffic and background sites for the period from 4 February to 7 March 2020. Both simulated and observed concentrations are presented as individual contributions from fossil fuel (FF), biomass burning (BB), and their combined total (total).

Figure 5 Hourly time series of observed (eBC) and simulated (BC) concentrations at urban traffic and background sites for the period from 26 December 2020 to 21 January 2021. Both simulated and observed concentrations are presented as individual contributions from fossil fuel (FF), biomass burning (BB), and their combined total (total). Please note the different scale for the y axis of the background and traffic stations.

This meteorological pattern contributed to the elevated eBC concentrations observed at both sites, particularly at the traffic station where levels reached up to 18 µg m⁻³ (Fig. 5). In these circumstances, the hybrid system struggled to reproduce the observed concentration pattern. A possible explanation is that it is difficult for the GRAL model to accurately simulate the PBL height over the urban area. As shown in Sect. 4.2.1, while GRAL generally agreed with the observations at 00:00 and 12:00 GMT, it failed to produce a realistic PBL height during certain sporadic episodes, particularly under neutral conditions at night. It can therefore be hypothesised that GRAL has difficulty simulating the very stagnant meteorological conditions that can occur in the Po Valley. Although this limitation exists, similar challenges have been documented in the literature for various other models and regions, as noted by Fay and Neunhäuserer (2006), Saide et al. (2011), Tominaga and Stathopoulos (2016), and Travis et al. (2022). Continuous measurements of the PBL height in Modena would be necessary to further analyse the accuracy of GRAL during the day with respect to this parameter.

Looking more closely at specific days, namely 3–4, 13–14, and 18–19 January 2021, significant peaks in BC levels were observed. At the traffic site, the peak concentrations on these days were 14.1, 18.5, and 16.3 µg m⁻³, respectively. Analysing the episodes on 13–14 and 18–19 January 2021, while direct BC data at the background site are not available for this period, the total PM₁₀ and PM_2.5 concentrations recorded at the same two locations (Fig. S7) confirm the presence of meteorological and emission conditions conducive to pollutant accumulation, with peak values of 105 µg m⁻³ for PM₁₀ and 66 µg m⁻³ for PM_2.5. Figure S2 shows a significant increase in atmospheric pressure at OSS, suggesting the establishment of a high-pressure system in the area. This atmospheric scenario was accompanied by calm and stable meteorological conditions, with an average wind speed of less than 2.5 m s⁻¹ (Fig. 3b). Furthermore, these conditions contributed to the formation of temperature inversions not only during the night but also during the day, as shown in Figs. S9 and S10. Taken together, these conditions probably facilitated the accumulation of pollutants near the ground, thereby limiting their vertical and horizontal dispersion. Comparable meteorological conditions were also observed on 16 January 2021. However, unlike the previously analysed days, 16 January fell on a Saturday, which is typically associated with lower anthropogenic emissions from traffic and industrial activities compared to weekdays. Consequently, the BC and PM concentrations on this particular day were lower than those recorded on 13–14 and 18–19 January 2021. On the other hand, the interpretation of the episode recorded on 3 and 4 January 2021 is somewhat more challenging. On these days, the total PM₁₀ and PM_2.5 concentrations were typical of the period, with no significant peaks (Fig. S7). In particular, daily mean PM₁₀ concentrations were around 30 and 25 µg m⁻³ at the traffic and background stations, respectively, while daily mean PM_2.5 concentrations were around 14 µg m⁻³ at the background station. Vertical temperature profiles derived from soundings at 00:00 and 12:00 GMT at SPC (Fig. S8) indicated only a shallow inversion layer near the ground, in contrast to the strong inversion layer observed on 13–14 and 18–19 January. In addition, other meteorological variables recorded at OSS did not suggest conditions conducive to a particular episode. It is possible that during the hours of peak concentration (18:00 GMT on 3 January and 08:00 GMT on 4 January), the inversion layer, combined with a temporary reduction in wind speed, may have decreased for a limited number of hours. This may have exacerbated the eBC concentrations recorded at high temporal frequency, while the daily averages of PM₁₀ and PM_2.5 were moderately affected by these short episodes. Furthermore, the presence of sporadic high-emission sources near the traffic station, such as high-emitting vehicles passing nearby or idling in nearby car parks, could have caused the peak on 4 January. Indeed, this peak was not observed at the background site and was dominated by fossil fuel contributions, suggesting the possible presence of local sources.

Despite the decrease in the performance of the hybrid modelling system compared to the first period, the metrics of FB, NMSE, FAC2, and NAD closely match the results reported by Lugon et al. (2021) when simulating BC concentrations in a suburban street network in Paris using the Street-in-Grid model. Furthermore, the corresponding metrics of FB, NMSE, FAC2, and NAD in this study for both periods, reported in Table 7 together with the previously mentioned indices, largely meet the acceptance criteria established by Hanna and Chang (2012) for urban dispersion modelling, where the benchmark values are defined as follows:

• FAC2 > 0.30

• NAD < 0.50

• |FB| < 0.67

• NMSE < 6.

Table 7 Statistical analysis of hourly total BC concentrations at urban traffic and urban background sites using the hybrid modelling system.

The MWAA model, applied to MA200 measurements with additional hypotheses about α for different sources, enabled us to discern the contributions of different sectors to the observed eBC concentrations (see Sect. 2.1). This approach allowed for the assessment of the hybrid system's performance in reproducing concentrations of two distinct sources: fossil fuel and biomass burning. Figures 6 and 7 present the diurnal patterns of observed and modelled concentrations, respectively, for the first and second simulated periods for total BC (eBC/BC; panel a), BC from fossil fuel (eBC-FF/BC-FF; panel b), and BC from biomass burning (eBC-BB/BC-BB; panel c). In addition, in order to assess the influence of emission uncertainties in the model outputs and explore potential outcome ranges, the figures for BC-FF and BC-BB also include the uncertainty ranges resulting from the emission scenarios described in Sect. 2.5. Error bars indicate the minimum and maximum average concentrations derived from different scenarios. The points correspond to the average concentrations of the base case scenario. Additionally, to provide additional details, blue-shaded areas depict the interquartile range derived from the base case scenario, while red-shaded areas represent the interquartile range of observed concentrations.

Due to a limitation in the MWAA model, it is not possible to separate the influence of non-exhaust and resuspension from the total eBC concentrations without additional experimental data. Therefore, where not explicitly stated, the modelled BC-FF concentrations in this section include the combined contribution of fossil fuels, traffic non-exhaust, and traffic resuspension.

During the first period (as shown in Fig. 6a), it is noteworthy that despite a slight underestimation of the measured concentrations occurring between 03:00 and 06:00 GMT and a slight overestimation observed between 18:00 and 21:00 GMT, both occurring at the urban background site, the simulation shows good agreement with measurements at both locations. Looking at the specifics of the eBC-FF/BC-FF concentrations (Fig. 6b), the base case scenario proves to be generally effective at representing the average trends, despite a slight underestimation between 18:00 and 21:00 GMT at both sites during the second peak of the daily rush hours. This discrepancy could probably be attributed to an inherent underestimation of traffic flows during this specific time period, a trend that is consistent with the results of the second period and with the findings of Veratti et al. (2021), who performed a comprehensive analysis of the same area with respect to nitrogen oxide concentrations. In order to test this hypothesis in depth, we compared the measured traffic flows from 45 induction loop sensors located at traffic lights in the vicinity of the two monitoring stations with the traffic flows simulated by PTV VISUM on the reference roads around these stations. The results showed that on 36 out of the 57 simulation days, the traffic flows between 18:00 and 21:00 GMT were underestimated, with NMBs of up to −18 % for roads near the urban traffic station and up to −13 % for roads surrounding the urban background site. Despite this underestimation, the overall agreement of model results with observations further confirms the effectiveness of a bottom-up approach in accurately representing traffic emissions.

Figure 6 Comparative analysis of modelled (BC, blue) and observed (eBC, red) mean daily cycles during the first period for total BC concentrations (a), fossil-fuel-derived BC (b), and biomass-burning-derived BC (c). The shaded regions in all three panels depict the interquartile range encompassing observations (red) and model outcomes (blue) for the base scenario. Error bars in panels (b) and (c) delineate the minimum and maximum average concentrations derived from the simulated scenarios in this study. Dots correspond to average concentrations from the base case scenario.

Focusing on the contribution of biomass burning to total BC (Fig. 6c), it is interesting to note that during the first half of the day, particularly between 03:00 and 11:00 GMT, the observations are close to the lower bound indicated by the error bars. This suggests that the collective average of all scenarios tends to slightly overestimate the observed trend during this period. Conversely, the experimental eBC-BB pattern could also be underestimated, possibly due to measurement uncertainties associated with the MA200 instruments, as reported in recent studies (Alas et al., 2020; Li et al., 2021). More specifically, the observations between 03:00 and 11:00 GMT show stronger agreement with the results generated by the DMH-2 scenario at the urban background location and by the DMH-4 scenario at the urban traffic location. Conversely, the peak in eBC-BB concentrations occurring after 19:00 GMT is represented remarkably well by the average concentration derived from all modelled scenarios. These results lead to two plausible interpretations. The first suggests that the hourly modulation profile used to distribute the emissions over the day takes into account a morning peak that may not be prevalent on average among the sources within Modena (see Fig. S12 for the emission daily modulation profile used in this study). The same hypothesis suggests that the distribution of emissions among the different districts of the city could possibly be an intermediate situation between those presented in this study (i.e. the DMH-1, DMH-2, DMH-3, DMH-4, and DMH-5 scenarios). The second alternative explanation suggests that the emission allocation is a mixture of those represented by the DMH-2 and DMH-4 scenarios and that the hourly modulation profile underestimates the emission peak during the night hours when the observations express the highest concentrations.

When analysing the second period, contrasting behaviour in the modelled concentrations is observed between the traffic and background sites (Fig. 7a). At the urban background station, the modelled concentrations are generally overestimated during most of the night and morning hours. However, the magnitude of the night peak is captured well, albeit with a 1 h delay, probably due to the emission modulation profiles for both FF and BB imposed on GRAL during the Christmas holidays (Fig. S12), which do not seem to fully represent the observed pattern. On the other hand, the measured concentrations at the traffic site are consistently underestimated during every hour of the day, although the hourly trend is represented well by the hybrid modelling system.

Figure 7 Comparative analysis of modelled (BC, blue) and observed (eBC, red) mean daily cycles during the second period for total BC concentrations (a), fossil-fuel-derived BC (b), and biomass-burning-derived BC (c). The shaded regions in all three panels depict the interquartile range encompassing observations (red) and model outcomes (blue). Error bars in panels (b) and (c) delineate the minimum and maximum average concentrations derived from the simulated scenarios in this study. Dots correspond to average concentrations from the base case scenario.

Despite the divergent results observed at the two stations, it is important to bear in mind that the concentrations represented at the two stations refer to two different periods: measurements were carried out from 26 December 2020 to 7 January 2021 at the urban background site and from 26 December 2020 to 21 January 2021 at the urban traffic site. Therefore, these datasets are not directly comparable (see previous explanations in this section). When analyses are restricted to periods common to both sites, model concentrations at the traffic site show similar behaviour those at to the urban background site. The hybrid modelling system reports an MB of 0.42 µg m⁻³, corresponding to +16 % of the NMB, while the Pearson correlation coefficient between modelled and observed values is 0.34, which is very similar to the statistics reported at the background site (Table 7). This overestimation is mainly due to three different episodes occurring between 27–28 December 2020, 4–5 January 2021, and 6 January 2021. The first and last episodes are evident at both sites, with both FF and BB affected by an overestimation of the same magnitude in relative terms. This suggests that the local meteorology was not captured well by GRAL, probably resulting in a lower PBL height than the actual values. On the other hand, the episode between 4 and 5 January 2021 is dominated by an overestimation of only the BB component at both sites. In this case, the concept of heating degree days used to assign the daily BB emissions to each simulated day may not have accurately represented the actual emissions for that specific day.

Looking more closely at the contributions of individual sources to the total $\mathrm{eBC}/\mathrm{BC}$ concentrations at the urban background site, a detailed examination shows that the overestimation in the morning peak is due to an overestimation of biomass burning emissions (Fig. 7c), which is consistent with the results of the first period. Conversely, the night peak is well captured by the average concentration derived from all biomass burning scenarios considered. The modelled BC-FF concentrations, although showing very limited variability throughout the day, are in good agreement with the measurements, suggesting that a precise modulation of traffic and other fossil-fuel-related sources for holidays (which predominantly characterise most days of the campaign at the background site) could have been applied.

At the traffic site, biomass burning concentrations are simulated well by the modelling system, with observations closely matching the averaged concentrations from all simulated scenarios for most hours except for the morning peak, confirming the previously developed hypothesis. In more detail, the results of the second period suggest a more plausible scenario where a single daily peak occurs at night between 20:00 and 22:00 GMT. Moreover, the daily modulation profile based on the concept of heating degree days seems to be quite effective at consistently reproducing biomass emissions for domestic heating at different winter temperatures, despite the overestimation reported between 4 and 5 January 2021. The modelled BC concentrations originating from fossil fuels at the traffic site show a tendency to converge towards the upper limits of the error bars corresponding to the upper range of the BC/PM speciation factors recommended by the EEA, i.e. by the combination of the TRF-exh-3, TRF-nexh-3, and TRF-bck-3 scenarios. While this result might suggest that the upper range of the speciation factors might reflect the composition of the vehicle fleet in Modena, we believe that the source of this underestimation is mainly related to the difficulties encountered by the modelling system in simulating specific meteorological events and the presence of emitting sources that pose challenges to the simulation.

As discussed earlier, instead of consistent underestimation of the observed trend throughout the simulated period due to low speciation factors, distinct episodes of significant underestimation can be identified (see Fig. 5). These episodes were characterised by robust meteorological stagnation (such as the events of 13–14 and 18–19 June 2021) that GRAL is unlikely to have captured. To support the hypothesis of inaccuracies in the simulation of the pronounced atmospheric stability observed on 13–14 and 18–19 January 2021, Fig. S11 in the Supplement provides an overview of the hourly PBL height simulated by GRAL in Modena, a critical variable for vertical mixing and pollutant dispersion. Although a direct comparison with measured values within Modena is not possible due to a lack of measurements, it is evident that the diurnal PBL height for the critical days mentioned (13–14 and 18–19 January 2021) is close to or even higher than that of the previous days, without any particular trend towards lower values. Consequently, the vertical dispersion of BC emitted at the city scale on these days was somewhat overestimated. Despite these inherent limitations, the results obtained from both the first and second periods show that the method used to estimate fossil fuel emissions, especially those from traffic, has proved its effectiveness in accurately reproducing the observed trend.

A limitation of our study is that the modelled BC concentrations are treated as inert and no external or internal mixing coatings were simulated by the hybrid modelling system. While the simulation of coating absorption is particularly challenging (Curci et al., 2019), especially for Lagrangian models (Doran et al., 2008; Holmes and Morawska, 2006), omitting this aspect could introduce potential uncertainties during model validation when assuming reference MAE values for aethalometer measurements. In fact, different mixing states can lead to an enhancement of the BC absorption efficiency (E_abs) compared to the theoretical pure form of BC. Studies carried out at different locations around the world have reported a wide range of E_abs values, from almost negligible (Cappa et al., 2012) to 2.4 (Peng et al., 2016), with most falling in the range of 1.2–1.6 (Moffet and Prather, 2009; Bond et al., 2006; Liu et al., 2017; Schwarz et al., 2008). To account for the uncertainty associated with MAE values when comparing modelled and measured concentrations, we investigated the impact of varying MAE. In particular, although we found that about half of the modelled BC concentrations in Modena originate from the city itself (as detailed in Sect. 4.2.3) where we assume negligible coating, especially when analysing concentrations close to emission sources, it is reasonable to expect that regional BC emissions originating from outside Modena may undergo ageing processes and thus experience absorption enhancement. Taking these factors into account, we examined the effects of adjusting the MAE within a range of ±20 %, which represents half of the average absorption enhancement range (1.2–1.6) documented in numerous studies in the literature. Figure S13 shows the linear regression between modelled and measured BC concentrations, using data from both simulated periods and the reference MAE values presented in Table S1. Setting the intercept to zero yields slope values of 0.83 and 0.65 when evaluating total BC concentrations at the background and traffic sites. At the same sites, the RMSEs of the regression model are 1.0 and 1.7, respectively. Focusing only on the contribution of biomass burning, the slopes of the regression lines are 0.98 and 0.77 at the background and traffic sites, respectively, with corresponding RMSE values of 0.7 and 1.0. For the fossil fuel contributions, the slopes of the regression lines are 0.49 and 0.56, and the RMSEs are 0.63 and 1.2 at the background and traffic sites, respectively. By increasing the MAE by 20 % compared to the values reported in Table S1, both the slope of the regression line and the RMSE of the linear model increase by about 17 %. Concurrently, the NMB between modelled and observed concentrations shifts from +7 % to −11 % at the urban background site and from −0.16 % to −30 % at the traffic site when assessing total BC concentrations. Conversely, a 20 % reduction in the MAE results in a reduction of approximately 25 % in both the slope of the regression line and the RMSE of the linear model. This variation in MAE results in an NMB of +34 % at the urban background site and +5 % at the traffic site when analysing total BC concentrations. The outcomes of the linear regression, together with the NMB values for the three scenarios analysed, are presented in Table S3. These results include assessments of total BC concentrations as well as contributions from biomass burning and fossil fuels.

4.2.3 Dispersion-modelling-based source apportionment

Once the robustness of the hybrid modelling system in reproducing the observed trends was verified, the second phase of the study employed the same system to identify potential sources influencing urban BC concentrations in Modena. The design of the modelling framework allows direct differentiation between contributions originating from within the city (labelled city in Figs. 8 and 9) and from sources located outside the urban area (labelled bck in Figs. 8 and 9). This distinction is achieved through adjustments made to the CHIMERE code, enabling real-time apportionment for non-reactive BC. In addition, the ability of the Lagrangian model GRAL to store concentration fields for selected sources allows direct apportionment of different urban emission sectors.

Figure 8 shows the average spatial contributions from each simulated source, which were calculated by analysing hourly maps from both the first and second simulated periods. In panel (a) of this figure, the upper row shows the average contributions of BB, FF, traffic non-exhaust, and resuspension emissions from urban sources (labelled city). Meanwhile, the lower row shows the average contributions from the same sources but originating from areas outside the region of interest and transported over longer distances (labelled bck).

Figure 8 Modelled breakdown of BC contributions, differentiating between fossil fuels (FF), biomass burning (BB), traffic non-exhaust emissions, and resuspension (a). The upper row represents the BC contributions from local sources within the city (city), while the lower row depicts contributions from sources located outside the urban area (bck). Panel (b) displays the total average BC concentrations, accounting for contributions from all these sources. Black rectangles in the graphs represent buildings. Please note the different scales used for panels (a) and (b).

When focusing on the contributions originating from within the city, these maps underscore the significant impact of BB on the overall levels of BC, as evident in the top-left graph of panel (a). Indeed, even though BC emissions are potentially absent or extremely low within the historical city centre, concentrations tend to spike in the surrounding area, reaching up to 2 µg m⁻³ from this contribution only. This can be ascribed to the increasing density and height of buildings as one moves from rural areas into the historical centre. These factors give rise to intricate urban canyons that foster atmospheric stagnation, thereby hindering the dispersion and dilution of pollutants. Regarding FF contributions from urban sources (centre top-left graph of panel a), it is evident that traffic is the primary source of FF combustion, as indicated by the concentration peaks of BC. These peaks exhibit noticeably high spatial gradients, particularly prominent in streets characterised by heavy traffic, such as the urban ring road encircling the urban area, the small section of highway in the bottom-left corner of the domain, and several other bustling urban streets. Furthermore, non-exhaust emissions and resuspension (the centre top-right and top-right graphs of panel a, respectively) also make a notable contribution to total BC concentrations, even if their absolute values are lower compared to those of BB and FF, especially in the vicinity of the aforementioned heavy-traffic streets.

In the lower row of panel (a), the contributions from long-range transport are displayed. Although these maps show more uniform concentration patterns when compared with those at the city scale, BB concentrations remain the dominant contributor among the four source categories considered here. This underscores that emissions outside the urban centre also exert a noteworthy influence on urban BC concentrations. On the other hand, Fig. 8b represents the sum of each single contribution from panel (a). This plot displays the ability of the hybrid modelling system to simulate urban concentration gradients originating from different source locations, and it also underscores the ability of the Lagrangian modelling suite to reconstruct the impact of urban canyons, which play a crucial role in the accumulation of pollutants in the most critical areas of the city centre.

Based on the results presented in Fig. 8 and considering a receptor location equivalent in scale to the GRAL simulation domain, the average concentration over all cells in this area was calculated for all emission sources analysed here. The reference value reported for each contribution represents the average of all simulated scenarios for the respective contribution (see Tables 2 and 3), and the associated uncertainty is the resulting standard deviation from these scenarios. A breakdown of the local urban contribution shows that BB accounts for 35 ± 15 %, FF for 9 ± 4 %, non-exhaust emissions for 4 ± 2 % and resuspension for 4 ± 3 %. It is also worth noting that FF is mainly attributed to traffic sources, with other FF sources accounting for less than 1 %. In contrast, according to the ARPAE local emission inventory, total BB emissions are roughly evenly distributed among open fireplaces (10 %), conventional stoves (8 %), high-efficiency stoves (8 %), and advanced stoves and boilers (9 %). Since there were no available data for a more detailed distribution of these emissions, we applied a consistent approach to spatially allocate all BB sources. As a result, the simulated BC concentrations align, with the same percentage breakdown as the emission inventory. These results indicate that the urban contribution accounts for 52 ± 16 % of the total average BC concentration, with the related uncertainty estimated using the root sum of squares method assuming that the uncertainties associated with each quantity are independent of each other, that the relationship between the quantities being summed is linear, and that the uncertainties follow a Gaussian distribution. Regarding long-range transport, it contributes 28 ± 1 % to BB emissions, 15 ± 2 % to FF emissions, 3 ± 2 % to non-exhaust emissions, and 2 ± 2 % to resuspension emissions, with uncertainties estimated similarly to the urban contributions. A summary of these percentages is depicted in the Sankey diagram presented in Fig. 9. Compared to the urban simulations, the uncertainty associated with long-range transport is lower, especially for BB. This reduction is due to the smaller discrepancies observed in the different simulated scenarios compared to the city-scale results. Specifically, when considering BB at a regional scale, changes in the spatial distribution of emissions yielded almost identical results to the assessment of concentrations in Modena. This suggests that when using a regional chemical transport model to simulate long-range transport to an urban area at a resolution of 3 km, the spatial distribution of emissions in the surrounding region becomes a less critical factor. Furthermore, the effect of minimum/maximum BC speciation factors seems to have a secondary effect on long-range transport.

Figure 9 Sankey diagram showing the percentage distribution of BC concentrations by source area and emission category.

4.2.4 Contribution of traffic sources

Building on the general source apportionment previously reported, the availability of detailed traffic information facilitated a more detailed investigation of the traffic contribution. As outlined in Sect. 2.5.1, five different sources were used to estimate traffic emissions. These sources include modelled traffic flows on a dense road network (Fig. 1), historical traffic counts from induction loop spires at key intersections, a comprehensive breakdown of the fleet composition, an assessment of annual kilometres travelled by each vehicle class, and speed-dependent EFs that take into account road and vehicle characteristics. Three types of traffic emissions were identified: exhaust, non-exhaust, and resuspension. Each type was simulated as a separate source in the GRAL model, allowing for a direct assessment of their individual contributions within the city. Additionally, dispersion simulations were conducted separately for the two representative vehicle categories (light and commercial vehicles) simulated by the traffic model. The results indicate that when the traffic contributions are aggregated, 50 % can be attributed to exhaust emissions, with light-duty vehicles accounting for the majority (38 %). The remaining 50 % is equally divided between non-exhaust and resuspension emissions, each contributing 25 % (Fig. 10).

Figure 10 Sankey diagram illustrating the percentage distribution of BC emissions from urban traffic sources.

Few studies in the literature have explicitly analysed the different contributions of exhaust, non-exhaust, and resuspension emissions for BC or EC, making a comprehensive comparison difficult. Among the studies that explicitly refer to BC, Lugon et al. (2021) investigated the dispersion of traffic emissions in a limited area in Paris using comparable emission tools and modelling techniques to those used in this paper, showing a very similar breakdown between the three components (approximately 50 % exhaust emissions and the remaining 50 % was the sum of non-exhaust and resuspension emissions). Other research has combined source apportionment methods with chemical analysis to distinguish the contributions of different sources. For example, Demir et al. (2022) investigated exhaust and non-exhaust traffic emissions in a road tunnel using EC fractions and the positive-matrix-factorisation method. Their results show that non-exhaust emissions contributed up to 37 % of the total carbonaceous species emitted by vehicles. In contrast, Karakaş et al. (2023) assumed that EC is only generated by fuel combustion or resuspended by traffic circulation, excluding contributions from tyre, brake, or surface wear. Therefore, the methodology and assumptions used to determine BC in non-exhaust emissions can significantly impact the final results, contributing to greater uncertainties.

A closer comparison of the present work with that of Lugon et al. (2021) reveals that the similarities between these two studies can be attributed to the EFs from the EMEP guidebook, which both studies used as a reference. However, it is important to critically note that these factors are primarily based on a literature review of publications from 1995 to 2002, which were largely based on indirect estimates of total suspended particulates (TSP) and PM₁₀. Recent reviews of direct measurements of tyre EFs have reported values approximately 4.5 times lower than the EMEP guidebook value for PM₁₀ and 42 times lower for PM_2.5 (Charbouillot et al., 2023; Harrison et al., 2021). Conversely, experimental laboratory tests conducted in the latest years (Kim et al., 2022; Lyu and Olofsson, 2020) confirm the reference BC EFs for brake wear reported in the EMEP guidebook, showing similar or even higher values with discrepancies of around 15 %.

For road wear EFs, the methods available in the literature are far from exhaustive and their estimates are of the lowest quality compared to tyre and brake wear emissions (Ntziachristos and Boulter, 2019). Estimates of this component have been provided by various authors (Woo et al., 2022; Amato, 2018; Gehrig et al., 2010; Thorpe and Harrison, 2008), but none of them explicitly refer to the BC or EC, making comparison with the EMEP guidelines difficult. Nevertheless, it should be noted that the contribution of road wear compared to the other sources is generally modest.

To explore the uncertainty associated with the attribution of the non-exhaust estimates reported in this study, an additional scenario was tested for each of the three components (tyre, brake, and road wear). For tyre wear emissions, a correction factor of 4.5 was applied to the TSP EFs, as suggested by Charbouillot et al. (2023), to align with recent direct measurements, resulting in an average BC EF of 0.38 mg km⁻¹ per vehicle. For brake emissions, the EFs were increased by 15 %, resulting in an average BC EF of 0.25 mg km⁻¹ per vehicle. As it was difficult to compare road wear emissions with the results of recent studies, the associated uncertainty was conservatively estimated by applying a correction factor of 4.5, similar to that used for tyre emissions. This resulted in BC EF values ranging from 0.35 to 0.02 mg km⁻¹ per vehicle. The results show that when these correction factors are applied, the contribution of tyre wear emissions decreased from 20 % to 5 %, while brake wear concentrations increased by 1 %. For road surface emissions, the contribution varies from 5 % to almost negligible when emissions are increased or reduced by 4.5, respectively.

The uncertainties reported in Fig. 10 for tyre, brake, and road wear are computed as the standard deviation of the reported scenarios, resulting in 11 % for tyre, 1 % for brake, and 2 % for road wear emissions. The overall uncertainty in the non-exhaust concentrations is then calculated using the root sum of squares method, as previously described. Similarly, the uncertainty attributed to the BC resuspension is represented by the standard deviation between the scenarios simulated following the methodology proposed by Amato et al. (2012) and those performed with the NORTRIP methodology (see Sect. 2.5.1).

In order to further explain the percentage breakdown of the subcategories in the exhaust emissions, the methodology used to calculate different contributions is extended. The representative EF for the two simulated vehicle classes (light and commercial vehicles) is calculated by averaging the EFs of all subcategories, weighted by a combination of the fleet composition and estimated annual mileage per vehicle class. This integrated approach allowed the composition of the registered fleet to be adjusted to reflect the actual average presence of vehicles on the road (see Fig. S14 for the fleet composition resulting from this calculation). The same approach is then used to calculate the contribution of each subcategory (vehicle type, fuel, and Euro standard) to the total emissions of either light- or heavy-duty vehicles. From this analysis, passenger cars (PC) are the main contributor to the exhaust component, accounting for 35 % of total transport emissions, while mopeds and motorcycles, light commercial vehicles, and buses and trucks contribute 3 %, 10 %, and 2 %, respectively. In addition, Euro 4 diesel passenger cars, following the aforementioned methodology, are responsible for the largest share of BC emissions in the urban environment, accounting for 19 % of total transport emissions. Although their EFs are lower than those of vehicles with older emission standards (Euro 3, Euro 2, etc.; see also Fig. S15), their share in the actual fleet composition (number of registered vehicles adjusted by annual mileage per vehicle class) makes this vehicle class account for a large contribution.

The uncertainty in the total exhaust contribution is estimated to be 16 % based on the standard deviation of the three simulated scenarios considering different BC speciation factors (see Sect. 2.5.1). On the other hand, Kouridis et al. (2010), suggested that the uncertainty in the exhaust EFs at a speed of 50 km h⁻¹ estimated following the tier 3 approach, is around 0.001 g km⁻¹ for both petrol and diesel passenger cars (corresponding to 53 % and 7 % of the weighted EFs used for the same categories in this study), around 0.06 g km⁻¹ (8 %) for commercial vehicles, and about 0.076 g km⁻¹ (48 %) for buses. No data are available for mopeds and motorcycles.

5 Conclusions

The present article discusses the outcomes of an integrated modelling and measurement approach designed to assess black carbon (BC) levels and sources in Modena, a city located in the Po Valley and chosen as representative of a typical medium-sized urban setting in the region. The eBC measurements were performed with two multi-wavelength micro-aethalometers of the Aethlabs MA200 series during two winter periods at two different urban sites: the first from 4 February to 7 March 2020 and the second from 26 December 2020 to 21 January 2021. One site was selected to represent typical urban background conditions, while the other was specifically chosen to capture high traffic intensity. In parallel to the measurement campaign, we implemented a hybrid modelling approach integrating the NINFA modelling suite (built upon the chemical transport model CHIMERE) with the Lagrangian model GRAMM-GRAL. This approach aimed to provide spatially resolved insights into BC concentrations and to identify potential sources contributing to BC levels in the urban environment of Modena. GRAMM-GRAL was used to reconstruct BC concentration fields at a high resolution (4 m), taking into account local sources and the presence of obstacles (e.g. buildings, bridges, or portals) in the flow field reconstruction. Conversely, NINFA was used to estimate the influence of external sources within the urban environment of Modena on BC concentrations.

Absorptions estimated with the MA200 instruments were further partitioned using the multi-wavelength absorption analyser (MWAA) model, allowing differentiation between different components, namely FF (emitted from fossil fuel sources), BB (from biomass burning), and BrC (brown carbon). The same partitioning was also used to differentiate eBC concentrations between fossil fuels (eBC-FF) and biomass combustion (eBC-BB). The absorption related to BrC showed a consistent and remarkably low presence throughout the day in both periods (around 0.5 Mm⁻¹ in the first, 1.2 Mm⁻¹ in the second). On the other hand, FF, a marker for traffic-related emissions, showed distinct diurnal peaks during rush hours, in the morning (07:00–08:00 GMT) and in the evening (18:00–20:00 GMT) for both periods, with values up to 187.2 Mm⁻¹ (corresponding to 18.5 µg m⁻³ of eBC-FF concentration) at the urban traffic site. In contrast, BB showed less diurnal variability, with a minimum in the early afternoon (15:00–17:00 GMT) and a peak in the evening (20:00–21:00 GMT), and generally lower absorption compared to FF. A comparison of the two seasons showed significant differences. The period from 26 December 2020 to 21 January 2021 had higher FF and BB concentrations, likely due to different meteorological conditions such as lower wind speeds and increased atmospheric stability, leading to frequent thermal inversions and a higher frequency of stable atmospheric conditions. In addition, lower temperatures likely led to more biomass burning for heating, which increased eBC concentrations.

In the second part of the study, the concentrations reproduced by the hybrid modelling system were compared with the observations recorded by the MA200 instruments. During the period from 4 February to 7 March 2020, the hybrid system showed robust performance, with an MB of 0.02 µg m⁻³ (+2 %) at the urban background site and −0.12 µg m⁻³ (−5 %) at the traffic site. The linear correlation coefficients were 0.51 and 0.62 respectively. Importantly, the daily concentration peaks were consistently captured under varying meteorological conditions, including periods of high wind speed. However, during the period from 26 December 2020 to 21 January 2021, contrasting behaviour was observed between the traffic and background stations. At the urban background site, measurements were only available for the period from 26 December 2020 to 7 January 2021, which included 7 holidays out of 13 measurement days. Under these conditions, the model showed a tendency to overestimate concentrations from biomass burning emissions, especially during the early part of the day. This effect was similar to that observed in the first period but more pronounced due to the increased influence of the holidays. For this station, the MB increased to 0.25 µg m⁻³ (+15 %) and the correlation coefficient decreased to 0.38. Conversely, measurements at the traffic station covered the entire period from 26 December 2020 to 21 January 2021. Throughout this period, the meteorological conditions were favourable for the accumulation of pollutants, characterised by the presence of high-pressure systems and thermal inversions. In such circumstances, the model struggled to fully capture the complex meteorological dynamics, leading to a notable reduction in the mean bias (−0.92 µg m⁻³; −22 %) and correlation coefficient (0.34). In addition, the presence of occasional high-emission sources in close proximity to the traffic monitoring site may have contributed to increased eBC concentrations. Despite the decrease in performance during the second period, the statistical metrics considered in this evaluation, including FAC2, NAD, FB, and NMSE, met the established acceptance criteria for urban dispersion modelling. Furthermore, these results are consistent with metrics from similar studies simulating BC concentrations in urban areas.

The same system was used to identify the impact of different sources on the spatial average BC concentrations, both inside and outside the urban area of Modena. The analysis showed that the city itself contributes 52 ± 16 % of the total average BC concentration, while background sources account for about 48 ± 4 % of the local concentrations. A detailed breakdown of the city's contribution showed that BB is responsible for 35 ± 15 %, FF for 9 ± 4 %, non-exhaust emissions for 4 ± 2 %, and resuspension for 4 ± 3 %. In particular, FF is mainly attributed to traffic sources, with other FF sources playing a negligible role. On the other hand, BB emissions are evenly distributed between open fireplaces, conventional stoves, high-efficiency stoves, and advanced stoves, as detailed in the ARPAE local emission inventory. As for long-distance transport, it contributes 28 ± 1 % to BB emissions, 15 ± 2 % to FF emissions, 3 ± 2 % to non-exhaust emissions, and 2 ± 2 % to resuspension emissions. Finally, when analysing the traffic-related concentrations, we found that 50 % can be attributed to exhaust emissions (in particular 19 % of the total to Euro 4 diesel passenger cars), with the remaining 50 % divided equally between non-exhaust and resuspension sources, each contributing 25 %. However, it is important to emphasise that significant uncertainties remain in the estimation of non-exhaust emissions, particularly for the fraction attributed to tyre wear. Different scenarios tested in this study show that the application of correction factors derived from the recent literature can significantly change the estimates. Specifically, the contribution of tyre wear varies from 20 % to 5 %, brake wear varies from 3 % to 4 %, and road surface emissions vary from 5 % to almost negligible.

In future developments, this research can serve as a valuable basis for more-comprehensive air quality management strategies in Modena and in similar medium-sized urban areas in the Po Valley. Ongoing monitoring and assessment of BC concentrations, together with a deeper understanding of the sources contributing to BC levels, can provide valuable insights for policy makers and urban planners seeking to reduce the impact of BC on air quality and public health in the region.

Appendix A: Symbols and acronyms

Description of symbols and acronyms used in the text

Symbol	Description
α	Absorption Ångström exponent
babs	Aerosol absorption coefficient
$b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{BB}}$	babs from BC by biomass burning
$b_{\mathrm{abs},\mathrm{BC}}^{\mathrm{FF}}$	babs from BC by fossil fuel combustion
$b_{\mathrm{abs}}^{\mathrm{BrC}}$	babs from brown carbon
BB	Biomass burning
BC	Black carbon
BrC	Brown carbon
CCN	Cloud condensation nuclei
Cref	Multi-scattering correction factor
CTMs	Chemical transport models
EC	Elemental carbon
EEA	European Environmental Agency
EF	Emission factor
FF	Fossil fuel
IN	Ice nuclei
IR	Infrared
MAAP	Multi-angle absorption photometer
MAE	Mass absorption efficiency
MWAA	Multi-wavelength absorption analyser
NMVOCs	Non-methane volatile organic compounds
PBL	Planetary boundary layer
PC	Passenger car
PM	Particulate matter
PTFE	Polytetrafluoroethylene
RB	Rural background
SNAP	Selected nomenclature for air pollution
UB	Urban background
UT	Urban traffic
UV	Ultraviolet
WHO	World Health Organization