Aerosol absorption using in situ filter-based photometers and ground-based sun photometry in the Po Valley urban atmosphere

Light-absorbing aerosols (LAAs) are short-lived climate forcers with a significant impact on Earth's radiative balance. LAAs include dust aerosols, black carbon (BC) and organic light-absorbing carbonaceous aerosol (collectively termed brown carbon, BrC), which have also been proven to be highly toxic. In this study, aerosol absorption at five wavelengths (ranging from ultraviolet to infrared) was monitored continuously using filter-based photometers during two winter seasons in 2020 and 2021 in the city of Modena (southern central Po Valley, northern Italy), at two regulatory air quality monitoring sites, along with other pollutants (coarse particulate matter, PM₁₀; fine particulate matter, PM_2.5; O₃; NO; NO₂; and C₆H₆) and the vehicular traffic rate. The aerosol optical depth (AOD) and other column aerosol optical properties were concurrently monitored at four wavelengths by an AErosol RObotic NETwork (AERONET) sun photometer under urban background conditions within Modena. In situ absorption levels were apportioned to both sources (fossil fuel and biomass burning) and species (BC and BrC), while columnar absorption was apportioned to BC, BrC and mineral dust. The combined analysis of the atmospheric aerosol and gas measurements and of the meteorological conditions (in situ and from the ERA5 reanalysis) identified the location of potential urban sources of BC and BrC, most likely related to traffic and biomass burning. In situ data show different diurnal/weekly patterns for BrC from biomass burning and BC from traffic, with minor differences between the background and the urban traffic conditions. AERONET version 3 absorption aerosol optical depth (AAOD) retrievals at four wavelengths allowed the estimation of the absorptive direct radiative effect due to LAAs over the same period under the reasonable assumption that the AOD signal is concentrated within the mixing layer. AERONET retrievals showed a modest correlation of columnar absorption with planetary boundary layer (PBL)-scaled in situ observations, although the correlation improved significantly during a desert dust transport event that affected both in situ aerosol and columnar absorption, particularly in the blue spectrum range. A low correlation occurred between the contribution of BrC to aerosol absorption for the in situ and the columnar observations, with the BrC contribution being generally larger for in situ observations. Finally, evidence of a highly layered atmosphere during the study period, featuring significant spatial mixing and modest vertical mixing, was shown by ERA5-based atmospheric temperature profiles and by the large correlation of concurrent AERONET AOD retrievals in Modena and in Ispra (on the northwestern side of the Po Valley, ca. 225 km from Modena).

1 Introduction

Light-absorbing aerosols (LAAs) include dust aerosols and soot-like, graphitic, elemental carbonaceous light-absorbing particles, the latter of which are qualitatively named black carbon (BC). A wide range of experimental techniques, which rely on the different properties of LAAs, are available for the experimental measurement of BC. In order to harmonize the terminology used reporting the concentration of this species, the scientific community recommends reporting BC observations based on light absorption as equivalent BC (eBC; Petzold et al., 2013). eBC aerosol particles have a fairly constant refractive index across the ultraviolet–infrared (UV–IR) range (Moosmüller et al., 2009). The eBC concentrations are converted into the light-absorbing carbon mass concentration using the mass-specific absorption cross section (MAC). Another type of LAA is brown carbon (BrC; Andreae and Gelencsér, 2006; Laskin et al., 2015): the fraction of light-absorbing organic aerosol with optical properties that differ from those of BC, due to their enhancement in absorption towards UV wavelengths.

LAAs are short-lived climate forcers (∼ 1-week atmospheric residence time; Forster et al., 2021) and significantly affect the Earth radiative balance (Bond et al., 2013; R. Wang et al., 2016). In terms of global impact, BC has been shown to have a positive direct radiative effect at the top of the atmosphere (TOA) in the range of 0.71–0.82 W m⁻² (Chung et al., 2012; Bond et al., 2013; Lin et al., 2014). Estimates of global direct effect are lower for BrC, in the range of 0.04–0.57 W m⁻², than for BC (Feng et al., 2013; Lin et al., 2014; Saleh et al., 2014; Jo et al., 2016; Brown et al., 2018; Zhang et al., 2020). BrC concentrations are very spatially variable, and concentrations depend on the study specifics. Due to aerosol–cloud interactions, the overall effective radiative forcing of LAAs (i.e., the difference in their radiative effect between the present day and preindustrial times; Heald et al., 2014) is 0.15±0.17 W m⁻² for BC (Thornhill et al., 2021; Forster et al., 2021), with the largest part of this uncertainty mainly arising from the indirect and semi-direct effects exerted by aerosol on cloud condensation nuclei, ice nuclei and on the atmospheric lapse rate, along with the aerosol mixing state (Twomey, 1974; Charlson et al., 1992; Bond et al., 2013; Rosenfeld et al., 2014; Takemura and Suzuki, 2019).

In addition to aerosol effects on climate, the scientific literature has documented adverse effects on human health, which significantly affects life expectancy (Loomis et al., 2013; Cohen et al., 2017; West et al., 2016). The toxicological effect of particulate matter (PM) is known to depend on the aerosol size distribution and chemical composition (Pöschl, 2005). BC is one of the components with a proven harmful effect on human health (Janssen et al., 2012), and both long-term and acute exposure to increased eBC concentrations have been shown to increase the mortality risk (Ostro et al., 2015; Yang et al., 2021). Recent studies have also shown that exposure to increased eBC concentrations is positively associated with various health issues, such as ischemic heart disease and myocardial infarction (Luben et al., 2017; Magalhaes et al., 2018; Kirrane et al., 2019). In addition, Regencia et al. (2021) observed that short-term cumulative exposure to traffic-related eBC concentrations could adversely affect blood pressure, resulting in cardiovascular diseases. BrC has also been shown to have detrimental health effects, enhanced because of its enrichment in organic compounds (Chowdhury et al., 2019; Offer et al., 2022), possibly related to aerosol aging (Li et al., 2022; Tuet et al., 2017; Weitekamp et al., 2020).

The compilation of reliable and accurate emission inventories for eBC is critical for the development of robust air quality control strategies and the mitigation of global warming. However, the large uncertainty associated with source emission factors, PM speciation and eBC definition makes the implementation of systematic and harmonized emission estimates a challenging task. Despite these limitations, most studies identify road transport as the largest eBC emission source in Europe (Wang, 2015), followed by biomass burning and industry (European Environment Agency, 2013), as more recently confirmed by the analysis of the eBC emission change in Europe due to COVID-19 lockdowns (Evangeliou et al., 2021). Similar to BC, BrC can be directly emitted into the atmosphere during the combustion of fossil fuels, although its major source is biomass burning. BrC can also originate from secondary reactions, e.g., through aging processes or via the photooxidation of biogenic or anthropogenic volatile organic compounds (VOCs) (Laskin et al., 2015).

Several approaches have been proposed in the literature to measure LAAs, including photothermal interferometry, photo-acoustic spectroscopy, and online or offline filter-based light attenuation methods (Lack et al., 2014). The difference in the BC reported by these techniques increases when significant amounts of secondary organics are present (Kalbermatter et al., 2022). Both the interferometric and acoustic approaches can be considered to be thermally based measurements, as they quantify the fraction of absorbed optical energy that is rapidly transferred into the surroundings under a controlled light source emission. The main advantage of these techniques is their direct measurement of the absorption of particles while suspended in air; however, they both suffer from technical and operational limitations. For example, the photo-acoustic technique is very sensitive to atmospheric conditions, such as relative humidity, temperature and pressure (Langridge et al., 2013), while photothermal interferometry is sensitive to mechanical vibration, although it has recently gained new attention (e.g., Visser et al., 2020; Drinovec et al., 2022). Filter-based measurements are very simple to operate but have the main disadvantages of filter-related artifacts, such as filter loading and multiple-scattering effects within filter fibers and between the collected particles and the filter fibers, possibly leading to systematic errors in the measurements. With the aim of overcoming these limitations, different technical and analytical corrections have been developed to correct for the non-idealities of filter-based measurements (e.g., Weingartner et al., 2003; Petzold et al., 2005; Virkkula et al., 2007; Collaud Coen et al., 2010; Hyvärinen et al., 2013; Drinovec et al., 2015; Li et al., 2020) for filter absorption photometers common in field experiments and in air quality monitoring networks. The Aethalometer (Magee Scientific Aerosol, Ljubljana, Slovenia) is a commonly used filter-based photometer designed to measure LAAs at multiple wavelengths and at high temporal resolution, generally at fixed monitoring sites. Lightweight portable micro-Aethalometer instruments, such as the AE51 or the MA200 series (AethLabs, San Francisco, USA), have recently been developed and successfully used in complex urban environments for pedestrian exposure assessments (Viana et al., 2015; Good et al., 2017; Boniardi et al., 2021), mobile observations (Grivas et al., 2019; Liu et al., 2019, 2021), and vertical profile investigations via drones and balloons (Ferrero et al., 2011, 2014; Pikridas et al., 2019; Kezoudi et al., 2021). Despite their limitations, multiwavelength aerosol absorption observations by filter photometers have proven suitable for the application of source and component apportionment models, such as the “Aethalometer model” (Sandradewi et al., 2008) to apportion BC between wood burning and fossil fuel combustion emissions or the Multi-Wavelength Absorption Analyzer (MWAA; Massabò et al., 2015; Bernardoni et al., 2017) algorithm, which enables the disentanglement of the BC and BrC components of LAAs as well as the determination of their radiative forcing (Ferrero et al., 2021a).

Surface in situ aerosol measurements can provide important information about aerosol characterization and concentration for the lowest tropospheric layer. However, estimating the vertical distribution of aerosol particles or their columnar load remains crucial to completely understand their impact on the climate system. In order to meet this need, a worldwide network of calibrated sun/sky photometers – the AErosol RObotic NETwork (AERONET; Holben et al., 1998) – was developed, with the goal of measuring aerosol optical columnar properties, e.g., aerosol optical depth (AOD) and column single-scattering albedo (SSA). Numerous studies have attempted to compare in situ observations with ground-based columnar aerosol optical properties and have provided different results depending on the atmospheric mixing state, the aerosol vertical profile and the local/regional pollution conditions. Several authors have used the ratio between the surface in situ aerosol mass concentration or aerosol absorption and the boundary layer height (i.e., they rescaled surface data over this atmospheric layer) and have shown how this ratio underestimated sun photometry observations of AOD or absorption AOD (AAOD), respectively (e.g., Bergin et al., 2000; Slater and Dibb, 2004; Aryal et al., 2014; Chauvigné et al., 2016; Chen et al., 2019). These findings were consistent across various types of locations (e.g., rural background, moderately polluted or marine) and highlighted that, in those settings, the main factors generally limiting the representativity of surface in situ measurements of the atmospheric column are the aerosol mixing within the boundary layer (BL) and the presence of aerosol above the BL, which can contribute significantly to the extinction and absorption in the column.

Datasets allowing a worldwide trend analysis of LAA levels remain limited (Laj et al., 2020); however, according to both in situ (Collaud Coen et al., 2020) and ground-based columnar (Li et al., 2014) observations, in the Northern Hemisphere, particularly in the US and Europe, the aerosol absorption coefficient (σ_ap) has decreased over the last decade(s). More specific to the region of interest for our study, the Po Valley is a European hotspot for atmospheric pollution situated in northern Italy. Previous work focused on the Po Basin observed a decrease in both columnar AOD and in situ aerosol scattering and absorption in Ispra, on the northwestern side of the Po Valley, in the early 2000s (Putaud et al., 2014). This drop was consistent with a significant valley-wide decrease in in situ PM₁₀ and PM_2.5 ground levels (Bigi and Ghermandi, 2014, 2016), thanks also to a drop in primary PM emissions from vehicular transport. Similarly, a drop of ∼ 4 % yr⁻¹ was recorded over the period from 1997 to 2016 for the elemental carbon content in fog samples at the rural background site of San Pietro Capofiume (Gilardoni et al., 2020b).

Significant aerosol sources other than traffic remain present in the valley, e.g., biomass burning for domestic heating for several compounds (including organic aerosols and BC) and farming for NH₃ (a major PM precursor). Their role in PM levels was highlighted by the small decrease in PM across the basin (Ciarelli et al., 2021; Putaud et al., 2021) and in particle count in Modena (Shen et al., 2021) during the 2020 COVID-19 pandemic lockdown. Some studies in the Po Valley have addressed temporal and vertical variability in σ_ap in Milan, the largest city in the basin (Ferrero et al., 2011, 2014; Vecchi et al., 2018). These authors found a decline in BC levels within the mixing layer, with higher BC levels observed at the ground (i.e., 50–100 m) and a marked drop (more than 50 %) above the mixing height, with BC contributing to ∼ 10 % (∼ 8 %) of the overall PM₁ extinction (mass) at a Milan urban surface background site in winter. Other studies in the Po Valley focused on the effect of the reduction in NO_x and NH₃ on PM_2.5 levels (Veratti et al., 2023) as well as on the impact of biomass burning on surface aerosols, particularly at the rural background site of San Pietro Capofiume and the urban site of Bologna (Gilardoni et al., 2016; Costabile et al., 2017; Paglione et al., 2020). These latter studies highlighted the large absorption Ångström exponents (AAEs; Moosmüller et al., 2009) for biomass burning organic aerosol, ranging from ∼ 3 to 5, mainly due to aged aerosols in the aqueous phase and related to an increase in the organic aerosol $/$ BC mass ratio. Previous investigations of the spatial variability in PM surface observations highlighted the impact of large urban areas on the aerosol load, particularly for PM₁₀, using cluster analysis (Bigi and Ghermandi, 2014, 2016). A Europe-wide assessment of urban air quality by Thunis et al. (2017), based on a simplified dispersion model, estimated a 57 % contribution of in-city emissions to urban PM_2.5 in Milan, making this city the one with the largest self-contribution to local PM_2.5 across the European Union. Similarly the spatial variability in columnar aerosol load observed by ground-based remote-sensing instruments between Ispra and the Adriatic Sea, east of the Po Basin, showed larger AOD and lower SSA at the Ispra site (Clerici and Mélin, 2008), confirming the impact of in-valley combustion emissions.

Relying on these previous findings, the current study provides additional knowledge on LAAs in the Po Valley by investigating the temporal, spatial and columnar variability in σ_ap in Modena, an urban area representative of several cities in the basin. The city of Modena is located in the southern central part of the Po Valley. The study period is winter 2020–2021, and the experimental dataset includes both in situ and ground-based columnar observations. Additionally, source apportionment of σ_ap using the in situ and the ground-based columnar observations in Modena are compared to investigate the impact of low-level emissions and long-range transport on the aerosol optical properties, along with the first estimation of the LAA heating rate (HR) and its diurnal trend in Modena. Finally, more insight into the spatial and temporal variability in the different absorbing components in the Po Valley are provided by a comparison between columnar optical properties in Modena and Ispra. Below, we first describe the measurements that we will use and then address these topics.

2 Measurement site and methods

Modena (44.6^∘ N, 10.9^∘ E; 32 m a.s.l., above sea level; ∼ 180 000 inhabitants) is located on the central southern side of the Po Valley, northern Italy, and is a basin surrounded by the Alps and Apennine mountains. The basin area is affected by recurrent atmospheric temperature inversions in winter and low-wind conditions, leading to a buildup of atmospheric pollutants. As a result, the Po Valley is one of the largest European regions exceeding the daily PM₁₀ limits set by the European regulation (EC 50/2008) and by the World Health Organization (WHO) guidelines (WHO, 2021). The city is situated in a region with flat topography, 13 km north of the foot of the closest Apennine hills and 96 km south of the foot of the Alps, i.e., it is on the southern side of a wide (∼ 110 km) valley.

The latest bottom-up regional emission inventory for the Modena municipal area (ARPAE, 2020), reference year 2017, identifies traffic and domestic heating as the main PM₁₀ sources, contributing 38 % and 58 % of the total emissions, respectively, although Modena also hosts a few districts for light manufacturing (Selected Nomenclature for sources of Air Pollution, SNAP 3 and 4), contributing 3 % of total PM₁₀ emissions (Fig. 1). More specific to nonindustrial combustion (SNAP 2), most of buildings use compressed natural gas for both heating and cooking; thus, 99.4 % of SNAP-2 PM₁₀ emissions are estimated to be produced by biomass combustion for domestic heating (ARPAE, 2020).

Figure 1 Measurement site locations and settings. Areas outlined in purple indicate the manufacturing districts (© OpenStreetMap contributors 2023; distributed under the Open Data Commons Open Database License, ODbL, v1.0.).

As is common to most urban areas in the basin, vehicular traffic is the main source of NO_x emissions (78 % of total NO_x emissions; ARPAE, 2020), with a significant impact on local air quality (Ghermandi et al., 2020; Veratti et al., 2020) and on population exposure (Veratti et al., 2021). Modena's setting is quite representative of several midsize urban areas across the Po Valley, particularly in terms of traffic and domestic emission sources and topography.

2.1 In situ surface measurements

Two MA200 micro-Aethalometer instruments were installed in Modena, sampling from the gently heated ($\sim \mathrm{30}\pm \mathrm{2}$ ^∘C) glassware manifold inlet lines already in use for reactive gas monitors at the two regulatory air quality monitoring sites in town: Giardini (EoI code IT0721A; 44.637^∘ N, 10.906^∘ E; 39 m a.s.l.) and Parco Ferrari (EoI code IT1771A; 44.652^∘ N, 10.907^∘ E; 30 m a.s.l.). These two sites are representative of urban traffic and urban background conditions, hereafter referred as UT and UB, respectively (see Fig. 1 for their locations). The UT site faces a major four-lane (two lanes in each direction) road, with estimated median daily traffic counts of ∼ 20 000 vehicles, whereas the UB is within Modena's largest urban park at a distance of ∼ 120 m from the nearest road. The inlet height at both sites is approximately 4 m above the ground. The inlet has no size cut, i.e., the instruments are sampling total suspended particles.

The MA200 instrument is a filter absorption photometer that measures at five wavelengths (λ= 375, 470, 528, 625 and 880 nm) using polytetrafluoroethylene (PTFE) filter tapes. The AAE for in situ observations (hereafter AAEⁱ) was computed by a fit to absorption at all five wavelengths. The instruments were used in dual-spot sampling mode (1.09 and 1.10 firmware versions were installed during the study), and a compensation algorithm similar to the one proposed by Drinovec et al. (2015) was applied by the internal firmware. This firmware uses a multiple-scattering correction coefficient C^ref=1.3, which was chosen by the manufacturer in order to mimic the response of the AE33 Aethalometer (AethLabs, personal communication, 2020).

Aerosol absorption monitoring at the UT site was performed between 19 January and 21 April in 2020 and between 9 November 2020 and 8 March in 2021 at a time resolution of 1 min. At the UB site, aerosol absorption was monitored between 4 February and 13 October in 2020 at a 1 min time resolution and between 13 December 2020 and 20 March in 2021 at a 5 min time resolution. In order to compensate for the occasionally low absorption readings at the latter site, the 1 min raw transmittance counts at the UB site were firstly aggregated to 5 min and then used to compute the corresponding σ_ap via a transcription of the dual-spot compensation algorithm in the R programming language, as described in Drinovec et al. (2015). The MA200 measurements were screened depending on the status reported by the instrument. Flow calibration was performed before each filter change. The flow was set to 100 mL min⁻¹ in winter and increased to 125 mL min⁻¹ in summer due to the lower atmospheric concentrations. In the present study, only measurements from the winter months were analyzed, i.e., December, January and February. Strict lockdown restrictions in northern Italy due to the COVID-19 pandemic lasted from 8 March until 4 May 2020; therefore, the winter data reported here are representative of a business-as-usual scenario, partly spanning across two winter seasons. Absorption data were averaged to 1 h prior to the analysis in order to match the time resolution of the other analyzed variables.

A comprehensive uncertainty analysis for absorption observations by the MA200 has not yet been fully performed by the scientific community. Li et al. (2021) looked at the multiple-scattering uncertainty of its PTFE filter and suggested that multiple-scattering artifacts might lead to an overestimation of the absorption by BrC, with the bias dependent on the absorbing strength of the compound, i.e., a behavior qualitatively similar to that of the AE33 (Yus-Díez et al., 2021). Alas et al. (2020) performed a large intercomparison involving these devices, testing several MA200s (single-spot mode, 10 and 60 s time resolution) and highlighted a low unit-to-unit variability (ca. 2 %) across all wavelengths and good agreement (R²>0.93) for loading-corrected eBC when compared to the AE33. In the current study, an 8 % uncertainty was attributed to MA200 absorption, based on the mean standard error of the slope of the linear regression between the eBC by the MA200 and the AE33 found by Alas et al. (2020).

Regulatory air quality data were also available at the two sites and include NO, NO₂, O₃ (UB site only) and C₆H₆ (UT site only) at an hourly timescale. PM₁₀ and PM_2.5 (at UB site only) were also available, although at a daily timescale. The daily median (10th and 90th quantiles) PM₁₀ concentration at the UB site over the period from 2017 to 2021 was 24 µg m⁻³ (13 µg m⁻³ and 57 µg m⁻³), while the same statistics for PM₁₀ were 27 µg m⁻³ (14 µg m⁻³ and 63 µg m⁻³) at the UT site . Consistently, over the same period, hourly NO₂ at the UB site showed lower levels than at the UT site, with the two locations having a median (10th and 90th quantiles) of 23 µg m⁻³ (6 µg m⁻³ and 50 µg m⁻³) and 35µg m⁻³ (14 µg m⁻³ and 66 µg m⁻³), respectively.

Direct traffic counts were also available for the urban area during the period of investigation. These data were collected by 400 induction loops for traffic light control within the urban and suburban street network. Continuous vehicle counts from the induction loops nearest to the UT and UB sites were aggregated into 1 h total traffic data; these hourly aggregates were primarily used to highlight variability in traffic patterns. The uncertainty in the count from these devices is approximately 10 % (Bellucci and Cipriani, 2010).

Meteorological variables were provided by the regional weather monitoring network station within the urban area of Modena and include wind speed and direction (WS and WD, respectively), atmospheric temperature (T), relative humidity (RH), downward global radiation (Q) and atmospheric pressure (p). The site is on the roof of the municipality offices at 40 m above the ground and is the highest weather station in the urban area with data available over the study period. These data provide indications of the wind conditions inside the urban canopy but may differ from wind conditions at the 4 m height of the MA200 measurements.

2.2 Mixing layer height

Hourly estimates of the mixing layer height (MLH) were provided by the ERA5 reanalysis (Hersbach et al., 2023). The ERA5 reanalysis, provided by the European Centre for Medium-Range Weather Forecasts (ECMWF), results from the data assimilation of global observations into the Integrated Forecast System (IFS), a global numerical weather prediction model, to produce a globally complete and consistent dataset of physical quantities, continuous in time and space. ERA5 provides hourly estimates for several geophysical quantities, including the MLH, at a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid resolution over the period from 1940 to today. For this study, we extracted the MLH at the ERA5 grid point with the following coordinates: 44.75^∘ N, 11.00^∘ W. Modena is about 14 km south of this grid point, i.e., it lays between two ERA5 grid points, as the size of the ERA5 cell is ∼ 20 km × ∼ 28 km in this region. The MLH was extracted at the aforementioned grid point because it is representative of a cell over an area with flat topography, i.e., very similar to the area of Modena. The MLH estimates from ERA5 are used in the analysis because no experimental estimates of the planetary boundary layer height were available in town, and the closest location with regular atmospheric sounding (at 12:00 and 00:00 UTC) is in San Pietro Capofiume, a rural background site surrounded by flat topography, 53 km east of Modena. The ERA5 estimate of the MLH in Europe was assessed to be underestimated by ∼ 54 m (∼ 19 m) on average (median) based on a comparison between ERA5 and daytime radiosoundings by Guo et al. (2021). This underestimate represents a lower-end estimate, as soundings in Europe are generally taken around 12:00 UTC, i.e., when the MLH is quite developed.

It is worth noting that deficiencies have been observed in various planetary boundary layer and surface parameterizations for conventional meteorological models, e.g., IFS, the Weather Research and Forecasting (WRF) model, and the Consortium for Small-scale Modeling (COSMO) model (Martilli et al., 2021; Maroneze et al., 2021; Lapo et al., 2019; Battisti et al., 2017), leading to a challenge with respect to the characterization of strong thermal inversions (Mahrt, 2014; Acevedo et al., 2019), like those occurring in the Po Valley. All of these meteorological models typically rely on a single stability parameter, e.g., the Richardson number or the turbulent kinetic energy, to automatically estimate the MLH based on fixed thresholds, regardless of the wide range of possible atmospheric conditions, as these models cannot apply ad hoc methods for each individual situation. More specific to this study, MLH estimates from the IFS are based on the bulk Richardson number (Vogelezang and Holtslag, 1996), regardless of the atmospheric stability conditions, and the MLH is defined as the lowest level at which the bulk Richardson number reaches the critical value of 0.25 (ECMWF, 2017).

2.3 Ground-based columnar measurements

Column-integrated measurements of optical properties of the Modena urban atmosphere were collected by a multi-channel Cimel CE318 sun/sky photometer installed on the roof of the “Enzo Ferrari” Department of Engineering at about 20 m above the ground. The instrument is part of NASA's AERONET network (Holben et al., 1998). This site, within the grounds of the University of Modena and Reggio Emilia and representative of residential background conditions, is on the southeastern edge of the urban settlement, whereas the UB and UT sites are on the west side of the town at a distance of 4km and 3.5kmfrom the sun/sky photometer location, respectively (Fig. 1).

In our analysis of the Cimel data, we considered both Level-2.0 and Level-1.5 version 3 almucantar retrievals at four wavelengths (λ= 440, 675, 870 and 1020 nm) (Sinyuk et al., 2020). Level-2.0 absorption data are more robust (e.g., Dubovik et al., 2000), but Level-1.5 data provide more matches with surface measurements, as discussed below. The almucantar retrievals provide several columnar properties, including absorption aerosol optical depth (AAOD), SSA, the depolarization ratio and the lidar ratio at the four wavelengths, as well as the particle volume size distribution and AOD apportioned to submicron and supermicron aerosols from which the fine-mode fraction (FMF) is calculated (O'Neill et al., 2003). AERONET retrievals also allowed the estimation of the scattering $\text{AOD}\left(\text{SAOD}\right({\mathit{\lambda }}_j)=\text{Total AOD}({\mathit{\lambda }}_j)-\text{AAOD}({\mathit{\lambda }}_j\left)\right)$ for each of the four wavelengths, the wavelength dependence of SAOD, i.e., the column scattering Ångström exponent (hereafter SAE^c), and the column absorption Ångström exponent (AAE^c).

There is much information in the literature about the uncertainty in AERONET products (e.g., Eck et al., 1999; Andrews et al., 2017; Sinyuk et al., 2020; Kayetha et al., 2022). A fixed AOD uncertainty set to 0.01 was used, following Eck et al. (1999). For the current study, an AOD-dependent uncertainty in SSA at 440, 675 and 875 nm for AERONET v3 retrievals was estimated based on the data from an urban site in Sinyuk et al. (2020), ranging from 0.017 at 440 nm when AOD₄₄₀ is 0.7 to 0.103 at 870 nm when AOD₄₄₀ is 0.03. The overall uncertainty for the analyzed aerosol parameter, e.g., AAOD, was estimated as a propagation of the uncertainties of AOD and SSA and ranged between 0.011 and 0.033.

In addition, the direct radiative effect (DRE) at the top of the atmosphere (TOA) and bottom of the atmosphere (BOA), retrieved by AERONET under clear-sky conditions in winter, were considered. As the atmospheric aerosol is characterized by a significant absorptive capacity, the difference between the DRE at TOA and BOA (hereafter ΔDRE_atm) represents the instantaneous radiative power density absorbed along the atmospheric column by the aerosol within that atmospheric layer (Chakrabarty et al., 2012; Kedia et al., 2010). ΔDRE_atm is expressed in watts per square meter (W m⁻²), which is the common metric used in the literature to quantify the integrated radiative power density absorbed by the aerosol in the atmosphere (Kedia et al., 2010; Das and Jayaraman, 2011; Bond et al., 2013; Heald et al., 2014). However, as demonstrated in Ferrero et al. (2014), a more useful parameter is the absorptive DRE (ADRE) of atmospheric aerosol, which can be computed simply by normalizing ΔDRE_atm by the atmospheric thickness Δz hosting most of the LAA; this thickness in the Po Valley can be generally assumed to correspond to the MLH. The ADRE represents the radiative power absorbed by the aerosol per unit volume of atmosphere (W m⁻³). The advantage of using ADRE in the Po Valley environment in wintertime is that, in this case, most of the AOD signal is built up within the mixing layer, as shown by both Ferrero et al. (2019), who found that up to 87 % of the AOD signal in Milan was generated within the mixing layer, 8 % was generated in the residual layer and 5 % was generated in the free troposphere, and by Barnaba et al. (2010), who found similar figures at the Ispra background site. This means that, if the thickness Δz is the MLH, the ADRE will refer to the MLH with an expected maximum overestimation of approximately 13 % (i.e., roughly the amount of the aerosol optical depth above the MLH). From the ADRE, the instantaneous heating rate (HR, K d⁻¹) can be computed as follows (Ferrero et al., 2014):

$$\begin{array}{}\text{(1)}& \text{HR}={\displaystyle \frac{\text{ADRE}}{\mathit{\rho }{C}_p}},\end{array}$$

where ρ is the air density and C_p (1005 J kg⁻¹ K⁻¹) is the isobaric specific heat of dry air. The most important advantages of this AERONET-based approach to derive the LAA HRs are as follows: (a) the possibility of obtaining a rapid HR estimation to investigate the HR temporal evolution during a selected time period and (b) the possibility of deriving the HR using a well-established network (AERONET), allowing a global comparison of the output. This approach is limited because HR can be obtained directly by the AERONET retrievals only if most of the AOD signal is built up within the mixing layer (thus with an expected overestimation of ∼ 13 %). Due to these limitations, the analysis of the HR is limited to retrievals collected on days without a significant dust content. This screening process followed the same process used for in situ data described in Sect. 2.4.1, as the apportionment of in situ data suffered from a similar limitation to that of the HR analysis.

Furthermore, in some of the analysis described below, the in situ and columnar data were compared, requiring temporal matching of the two datasets. An in situ–columnar observation match is considered successful when the AERONET retrieval occurred during the hourly averaged in situ measurement. Level-1.5 version 2 AERONET data are known to have large uncertainty when AOD at 440 nm (AOD₄₄₀) is less than 0.4 (Dubovik et al., 2000). In order to maximize the availability of columnar measurements for the analysis, Level-1.5 data were used. Level-1.5 data points with AOD₄₄₀≤0.2 were discarded from the analysis, and the data remaining after the AOD screening are referred to as “L1.5^*” in what follows.

2.4 Source apportionment of in situ and columnar data

Both in situ and columnar data were apportioned according to the aerosol spectral properties, i.e., solving the balance of the aerosol absorption based on its dependence on the AAE, on the absorbing species and on the wavelengths. Two different apportionment approaches were used for the in situ and the columnar observations, although based on the same foundation. The approach applied to the in situ data requires at least five wavelengths to ensure stability (Bernardoni et al., 2017).

2.4.1 In situ apportionment

In situ aerosol absorption σ_ap was apportioned to species (black carbon and brown carbon, referred to as ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC}}\left(\mathit{\lambda }\right)$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}\left(\mathit{\lambda }\right)$, respectively) and sources (fossil fuel and biomass burning combustion, referred to ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{ff}}\left(\mathit{\lambda }\right)$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BB}}\left(\mathit{\lambda }\right)$, respectively) using the Multi-Wavelength Absorption Analyzer (MWAA) model (Massabò et al., 2015; Bernardoni et al., 2017). This model assumes an equivalence between the absorption Ångström exponent (AAE; Moosmüller et al., 2009) of BC and that of fossil fuel (AAE${}_{\mathrm{BC}}^i=$ AAE${}_{\mathrm{ff}}^i$), and it assumes biomass burning to be the only source of BrC. Under these hypotheses, the MWAA model assumes that both the following equations hold for the total σ_ap(λ) at each wavelength:

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{ap}}={\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC}}\left(\mathit{\lambda }\right)+{\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}\left(\mathit{\lambda }\right)=\text{A}\cdot {\mathit{\lambda }}^{-{\mathrm{AAE}}_{\mathrm{BC}}^{\mathrm{i}}}+\text{B}\cdot {\mathit{\lambda }}^{-{\mathrm{AAE}}_{\mathrm{BrC}}^{\mathrm{i}}},\end{array}$$

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{ap}}={\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{ff}}\left(\mathit{\lambda }\right)+{\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{bb}}\left(\mathit{\lambda }\right)={\text{A}}^{\prime }\cdot {\mathit{\lambda }}^{-{\mathrm{AAE}}_{\mathrm{ff}}^{\mathrm{i}}}+{\text{B}}^{\prime }\cdot {\mathit{\lambda }}^{-{\mathrm{AAE}}_{\mathrm{bb}}^{\mathrm{i}}}.\end{array}$$

In Eqs. (2) and (3), AAE${}_{\mathrm{BC}}^{\mathrm{i}}=$ AAE${}_{\mathrm{ff}}^{\mathrm{i}}=\mathrm{1}$ was set, based on the AAEⁱ computed over five wavelengths during morning rush hour on winter weekdays at the UT site, consistent with fresh, uncoated BC particles (e.g., Liu et al., 2018). The AAEⁱ for BrC was determined by a preliminary nonlinear fit of Eq. (2), performed treating AAE${}_{\mathrm{BrC}}^{\mathrm{i}}$ as a free parameter (and resulting in an average AAE${}_{\mathrm{BrC}}^{\mathrm{i}}=\mathrm{3.9}$); AAE${}_{\mathrm{bb}}^{\mathrm{i}}=\mathrm{2}$ was set based on data from the literature for the Po Valley (Bernardoni et al., 2011, 2013; Costabile et al., 2017; Vecchi et al., 2018). A and B were then obtained for each sample via a multiwavelength fit of Eq. (2) (after fixing AAE${}_{\mathrm{BrC}}^{\mathrm{i}}$), and A^′ and B^′ were obtained via a multiwavelength fit of Eq. (3). It is noteworthy that the MWAA model neglects possible contributions from mineral dust. To limit uncertainties resulting from this, the days with a significant dust load were discarded prior the application of the MWAA model to the in situ data (i.e., whenever the in situ apportionment data are presented throughout the text, they have been screened for dust). Days with a significant dust content were first identified for the atmospheric column, using the particle volume size distribution estimated by the AERONET inversion (Sinyuk et al., 2020); the identification of dust events was performed qualitatively, based on the retrievals having a dominant coarse mode (e.g., Fig. S7b in the Supplement). These retrievals were subsequently double-checked by 72 h Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) back trajectories using Global Data Assimilation System (GDAS) 1^∘ resolution wind fields. Additionally, the impact of dust at the ground level was assessed based on the daily PM_2.5 / PM₁₀ ratio from the in situ measurements (Fig. S1), with a ratio of  ≤0.5 as a qualitative threshold for a dust event. For reference, the daily PM_2.5 / PM₁₀ ratio in winter between 2017 and 2021 at the UB site had a median of 0.71 and a 10th (25th) quantile of 0.53 (0.62), i.e., the two aerosol fractions are quite similar, as previously observed at most UB sites across the basin (Bigi and Ghermandi, 2016).

2.4.2 Columnar apportionment

The AAOD was apportioned to BC, BrC and mineral dust using the approach proposed in Bahadur et al. (2012), i.e., by directly solving the system of Ångström equations (see Appendix A) using the AERONET almucantar L1.5^* retrievals. The system includes Eq. (A1), describing the additive contribution of the AAOD from each species to the total AAOD, and Eq. (A2), describing the exponential dependence of the AAOD on the wavelength. This apportionment method neglects the mixing state of absorbing species (i.e., the aerosol is assumed to be externally mixed) and assumes that the observed AAOD is representative of a well-mixed sample of these species. Bahadur et al. (2012) estimated globally valid ranges of the AAE^c and SAE^c for BC, BrC and dust, parameters needed to solve the system of AAE equations, based on long-term, worldwide AERONET observations (Level 2.0 version 2).

For the current study, a tailored estimate of the AAE^c values for Modena was performed, based on the full time series of AERONET retrievals in Modena (from January 2000 to June 2021). The classification of aerosol species (BC, BrC and dust), in order to estimate their AAE^c values, was performed by combining the approaches of Bahadur et al. (2012), Cazorla et al. (2013) and Shin et al. (2019). Cazorla et al. (2013) suggests threshold values for SAE^c and AAE^c across the 440–675 nm range (hereafter SAE1^c and AAE1^c, respectively), which were applied for a preliminary classification (Fig. S2). Shin et al. (2019) combined the particle linear depolarization ratio and the lidar ratio at 1020 nm into a dust ratio coefficient χ_d,λ, estimating the contribution of dust and non-dust aerosol to the AOD. Following Bahadur et al. (2012), in order to disentangle the spectral properties of fossil fuel and biomass burning aerosol, first AAE^c for dust was assessed using the full L1.5^* time series (259 data points), based on the conditions SAE1^c<1, AAE1^c>1.5 and ${\mathit{\chi }}_{\mathrm{d},\mathrm{1020}\mathrm{nm}}>\mathrm{0.8}$. As the major source of biomass burning in the Po Valley is domestic heating during winter, AAE^c for BC was estimated based on the full time series of summer L1.5^* retrievals (1752 data points). The conditions applied in this case were SAE1^c>1.2 and AAE2${}^{\mathrm{c}}/$ AAE1^c>0.8, with an index of 1 indicating the range 440–675 nm and an index of 2 indicating the range 675–880 nm (Bahadur et al., 2012). Then, the AAE^c for BrC was computed by solving the AAE equation system on the L1.5^* non-dust winter retrievals over the period from 2015 to 2022 (89 data points). The AAE^c values for BC and BrC are based on datasets with different sizes, as fewer retrievals are available in winter due to the shorter daytime duration and clouds; to limit the possible bias induced by this difference in sample size and by the potential presence of outliers, the median ± median absolute deviance values of the AAE^c for dust, BC and BrC were computed for both Modena and Ispra (see below) and are reported in Table 1. Table 1 also includes values of column AAE^c from the literature for different absorbing aerosol types for comparison.

Bahadur et al. (2012)Bahadur et al. (2012)Dubovik et al. (2002)Giles et al. (2012)Russell et al. (2010)Zhang et al. (2022)Mallet et al. (2013)Kayetha et al. (2022)Zhu et al. (2021)

Table 1 Summary of the columnar absorption Ångström exponent (AAE^c) for BC, BrC and dust from this work and other studies in the literature. Rows are organized by wavelength.

^a These values generically refer to “urban/industrial/polluted aerosol”, i.e., potentially from a mixture of BC and BrC.
^b GRASP: Generalized Retrieval of Aerosol and Surface Properties (Dubovik et al., 2014).
^c These values generically refer to “dusty sites”.
^d SKYNET is a worldwide network of sun/sky photometers (Nakajima et al., 2020).

To assess the representativity of the AAE^c values derived for Modena, sun/sky photometer retrievals in Modena were also compared to AERONET data from Ispra (45.80^∘ N, 8.63^∘ E; 220 m a.s.l.; 225 km northwest of Modena) collected by a second Cimel CE318 sun/sky photometer within the AERONET network. Ispra exhibited a SAE${}^{\mathrm{c}}/$ AAE^c matrix very similar to that observed in Modena (not shown). The resulting AAE^c values for BC, BrC and dust in Modena and Ispra (Table 1) are consistent with most of the existing literature and the variability reported therein (e.g., Russell et al., 2010; Bahadur et al., 2012; Giles et al., 2012; Kayetha et al., 2022).

Finally, with reasonable confidence in the tailored AAE values for the different absorbing components, each AERONET retrieval at Modena and Ispra was apportioned by summarizing the solutions of the equation system as described by Bahadur et al. (2012). The apportionment was performed using the following two-step procedure, based on the assumption that AAE^c followed a normal distribution featured by the parameters in Table 1:

• step 1 – random extraction of AAE^c for all species at all wavelengths;

• step 2 – direct solution of the system of Ångström equations.

Steps 1 and 2 were repeated 10⁴ times for each retrieval in order to develop statistics of the AAE^c combination that provides a solution to the system. The time series of median AAE^c values was fairly stable over the measurement period at both sites, except during an intense episode of dust transport during which the AAE^c for BrC increased significantly and AAE2^c for dust dropped (Fig. S3). Both AAE^c for BrC and AAE2^c for dust values were on the tails of their respective distributions. It is worth noting that AAE^c refers to BC and not to eBC, as it results from a direct estimate of the absorption wavelength dependence of aerosol particles while suspended in air.

3 Results

3.1 Diurnal patterns for the in situ data

Figure 2 shows the medians and interquartile ranges of atmospheric species obtained from in situ observations along with the hourly traffic count from the induction loops closest to each monitoring site for winter (December, January and February) from early 2020 until March 2021. These data are screened for days with a non-negligible dust load, as specified in Sect. 2.4.1. The σ_ap at 528 nm for winter weekdays (Monday through Friday) and winter holidays (i.e., Sundays and local and national holidays) is in the top panel of Fig. 2 and represents the absorption by aerosol at about 4 m above the ground. Saturdays are excluded due to their mixed signal between a holiday and a weekday. The pattern of absorption apportionment components at 880 nm is also shown, followed by NO, NO₂, O₃ (UB only) and C₆H₆ (UT only). Figure 3, based on the same dataset as Fig. 2, displays the share of σ_ap at 375 nm due to BC from fossil fuel, due to BC from biomass burning and due to BrC, along with the variability in AAE over the 375–880 nm range. The medians and interquartile ranges for meteorological variables over the same period are shown in Fig. 4, while the hourly wind rose is shown in Fig. S4. Overall median and interquartile ranges for the dataset in Fig. 2 (i.e., with dust screening) and Fig. 4 are shown in Tables 2 and S1 in the Supplement, respectively.

Figure 2 Diurnal pattern of the median and the interquartile range values for total σ_ap at 528 nm, the apportioned σ_ap at 880 and regulatory gas compounds at the urban background (left) and urban traffic (right) air quality monitoring sites on weekdays (orange) and holidays (gray).

Figure 3 Diurnal pattern of the median and the interquartile range values for apportioned σ_ap at 375 nm and of the AAE (375–880 nm) at the urban background (a, c, e, g) and urban traffic (b, d, f, h) air quality monitoring sites on weekdays (orange) and holidays (gray).

Figure 4 Diurnal pattern of the median and the interquartile range values of meteorological variables in winter (December–January–February) during weekdays and holidays. All variables stem from an urban meteorological station except the mixing layer height, which is provided by the ERA5 reanalysis.

Table 2 Summary of atmospheric species and the traffic volume for the urban background (UB) and the urban traffic (UT) sites, showing the median (Med) and the 25th and 75th quantile (25Q and 75Q, respectively) values.

Table 3 reports a comparison of aerosol and BrC absorption reported by this and other studies at a few Po Valley sites. Values in Modena are generally higher than those reported in urban background Milan and the rural background sites of Motta Visconti and Ispra (Ferrero et al., 2021b; Gilardoni et al., 2020a; Laj et al., 2020; Zanatta et al., 2016). The data from these earlier studies come from filter absorption photometers, either Multi-Angle Absorption Photometer (MAAP) or Aethalometer instruments, with the latter instrument being corrected for multiple-scattering-induced bias based on co-located observations. For example, in Gilardoni et al. (2020a), a C_ref=3.0 based on Collaud Coen et al. (2010) was used to correct AE22 absorption. Compared with other southern European urban sites, Modena recorded larger σ_ap at 660 nm than Barcelona (Ealo et al., 2018) but lower σ_ap at 375 nm than Athens (Greece), mainly due to the large impact of biomass burning emissions in this city (Liakakou et al., 2020; Katsanos et al., 2019). A similar pattern is observed for BrC absorption. No MAAP co-location was available for the two MA200 instruments in Modena, leading to a larger uncertainty in their absolute readings; however, these two units showed good agreement with a MAAP during a BC intercomparison in urban background Athens (Stavroulas et al., 2022), where they exhibited a linear slope of 1.00 (R²=0.92) in winter and 1.07 (R²=0.92) in summer.

Gilardoni et al. (2020a)Gilardoni et al. (2020a)Ferrero et al. (2021b)Gilardoni et al. (2020a)Gilardoni et al. (2020a)Kaskaoutis et al. (2021)Zanatta et al. (2016)Ealo et al. (2018)Laj et al. (2020)Gilardoni et al. (2020a)Gilardoni et al. (2020a)Kaskaoutis et al. (2021)

Table 3 Summary of the mean ± standard deviation of absorption for aerosol and BrC based on this work and other studies in the literature. Rows are organized by wavelength. UB, UT and RB stand for urban background, urban traffic and rural background, respectively.

^a The 95 % confidence interval of the mean.
^b The geometric mean ± geometric standard deviation.
^c The median (10th–90th quantile).
^d The BrC determined on methanol extraction.

The diurnal pattern in the absolute levels of fossil fuel combustion species (Fig. 2) exhibits a similar pattern. There is an initial increase during the morning rush hour (08:00–10:00 LT, local time), followed by a drop at midday due to the dilution induced by an increased MLH depth despite the steady traffic rate. A second increase occurs at 18:00 LT, followed by a drop at approximately 20:00 LT that is delayed compared with the drop in traffic, possibly because of the shallow MLH in the evening or because of frequent thermal inversions at ground level. More specifically, on weekday evenings, ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ peaks at 20:00 LT, 1 h later than on holidays, at both UB and UT, with the former site recording ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ levels higher in the evening than in the morning.

This main pattern is followed by all atmospheric species except for secondary pollutants (e.g., O₃) and aerosols related to biomass burning (e.g., ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{bb}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$). The pattern features higher concentrations on weekdays than during holidays at both sites, as shown by the 7 %–240 % increase on weekdays, similar to the 230 %–250 % increase in traffic, confirming a major and local fossil fuel combustion source of direct origin. For most of these species, the absolute interquartile range (IQR) is larger on weekdays than during holidays. For gas-phase compounds, the IQR increased on weekdays, from 8 % for NO₂ at the UT site to 160 % for NO at the UB site. For absorption, the largest increase in the IQR (a 56 % increase) occurred for ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ at the UT site. The difference in variability between weekdays and holidays might be partly driven by the larger variability in fossil fuel combustion emissions during the former, along with the larger count of weekdays compared with holidays in the statistics.

${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{bb}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ exhibit a diurnal pattern featuring an increase from 17:00 to 23:00 LT, possibly triggered by the decrease in both the MLH and T, leading to the condensation of semivolatile organics and to an increase in biomass burning emissions. A similar diurnal pattern for BrC absorption was observed in UB Milan and RB Po Valley sites (Gilardoni et al., 2020a) as well as in the UB in Athens (Liakakou et al., 2020; Kaskaoutis et al., 2021). The weekly pattern for these two species is larger at the UT site, with an increase in the overall median values of ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{bb}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ of 22 % and 35 % during holidays, respectively, along with an increase in their respective IQRs of 16 % and 28 %.

The share of absorption at 375 nm for the three apportioned species (Fig. 3) shows a distinct diurnal and weekly pattern at the UT site, with ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{bb}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ being ca. 37 % larger during holidays, in contrast to ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ which is 32 % larger on weekdays. The UB exhibited a similar pattern, although with lower intensity. The holiday increase in biomass burning aerosol is probably linked to the longer period spent at home compared with weekdays and to a large recreational use of biomass burning in town, where most houses use compressed natural gas for domestic heating and cooking (99.4 % of SNAP-2 PM₁₀ emissions in Modena are from biomass combustion for domestic heating according to ARPAE, 2020). The diurnal pattern of the share of absorption from ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ at 375 nm is similar to the diurnal traffic count cycle, exhibiting larger values during weekdays and at the UT site. This supports the results of the apportionment and the hypothesis of the role of the MLH in the evening enhancement of absorption.

O₃ exhibits a “weekend effect” (Cleveland et al., 1974), common to most urban areas in Europe with a VOC-limited regime, i.e., ozone rises earlier in the morning on holidays due to the lower NO_x levels, leading to a more efficient photocatalytic cycle, and drops later in the evening due to the (later) increase in NO_x.

Atmospheric heating by aerosols based on σ_ap values in Modena was estimated by determining the HR from AERONET data, as detailed in Sect. 2.2. Figure 5a shows the complete HR time series obtained over Modena during the investigated period. Under an average (standard deviation) irradiance value of 386 (143) W m⁻², the average (standard deviation) HR was 1.61 (1.58) K d⁻¹. This value is consistent with data from Milan for wintertime, under clear-sky conditions, where a mean (± mean confidence interval) of 1.68±0.04 K d⁻¹ was found, when the incoming radiation was similar (441±148 W m⁻²). This latter point is important, as AERONET data are mainly available under clear-sky conditions; thus, the obtained HR data represent the upper limit for the site. In the Po Valley, the HR was shown to decrease by ∼ 12 % for every okta of sky covered by clouds (Ferrero et al., 2021b). With respect to the HR diurnal pattern, Fig. 5b shows the mean diurnal pattern of irradiance and HR under clear-sky and cloudy conditions. The incoming radiation peaked at 529±55 W m⁻² at 13:00 LT (Fig. 5b), while σ_ap peaked between 08:00 and 10:00 LT (Fig. 2). This causes an asymmetric HR diurnal pattern, characterized by a fast increase to the maximum at 11:00 LT (1.83±0.84 K d⁻¹) and a subsequent slower decrease till sunset (Fig. 5b), as is common under clear-sky conditions (Ferrero et al., 2018, 2021b).

Figure 5 Time series (a) and diurnal pattern (b) of the aerosol heating rate (HR) at Modena from AERONET retrievals.

The wind pattern in Modena features a mild mountain–valley breeze system along the Po Valley longitudinal axis, superimposed on the local wind circulation. During the investigated winter period, calm wind conditions (speed lower than 1 m s⁻¹) occurred 25 % of the time and an overall wind speed average of ca. 1.5 m s⁻¹ was recorded. Northwesterly winds, blowing from the higher side of the valley, dominated during daytime hours (11:00–17:00 LT) and were associated with the highest wind speed (occasionally above 9 m s⁻¹). The rest of the day featured local west–southwesterly low winds (wind speeds lower than 3 m s⁻¹) and some easterly winds (Fig. S4).

A conditional bivariate polar function was applied to NO_x, ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ at 880 nm and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ at 375 nm at both the UT and UB sites to identify the position of potential emission sources (Fig. 6) on weekdays and holidays, excluding days with a significant dust load at the ground. Wind speed and direction data were combined with atmospheric compounds' levels by the use of conditional bivariate polar functions (Uria-Tellaetxe and Carslaw, 2014), as implemented in the openair package for the R software (Carslaw and Ropkins, 2012). This tool provides information on both the direction and the distance of the (relatively local) emission sources that are contributing significantly to the observed concentration levels as well as on the wind direction sectors that might provide clean air masses.

Figure 6 Bivariate polar function applied to NO_x (a), ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ at 880 nm (b) and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ at 375 nm (c) at the UB (left panels) and UT (right panels) sites for weekdays and holidays.

At the UB site, NO_x and absorption exhibit slightly different directional patterns: both show an increase associated with slow south–southwesterly winds, particularly on weekdays, while ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ exhibits an increase on weekdays during northeasterly moderate winds, probably linked to traffic on the busy road 400 m in that direction. ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ exhibits larger values during holidays during southerly winds, which is different from the pattern for BC from fossil fuels and NO_x for the same period. At the UT site, the directional pattern between NO_x and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ is quite similar, highlighting the role of nearby traffic during weekdays and of the major east–west road south of the UT site which contributes mainly during holidays. Moreover, at the UT site, ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$ is higher during holidays and under southerly wind conditions. This latter increase occurs during evening/night hours (not shown), consistent with biomass burning from domestic heating for recreational use, with the increase probably enhanced by nighttime atmospheric stagnation. Finally, northwesterly moderate winds are associated with low levels of NO_x, ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$, mainly because northwesterly winds occur primarily at midday during maximum atmospheric mixing.

3.2 Comparing absorption optical depth from remote-sensing and in situ values

In situ and columnar aerosol optical properties were compared to assess (1) how representative the surface in situ aerosol optical measurements are of the mixed layer and (2) how the absorption within the MLH compares to the atmospheric column. Urban in situ and column data were compared over the whole time period, although simultaneous observations were mainly only available in February 2021. For this comparison, the in situ σ_ap values were rescaled over (i.e., multiplied by) the MLH height; this estimate of the integral aerosol absorption height corresponds to the case of vertically homogeneous σ_ap across the MLH. The in situ σ_ap in the IR spectral range (λ=880 nm) rescaled over the MLH height was generally larger than the AAOD (Fig. 7), for both the L1.5^* and L2.0 AERONET inversions, with mean normalized errors of MNE =2.2 and MNE =1.7, respectively. A better agreement occurred for blue wavelengths (λ= 470 and 440 nm for the in situ and the columnar observations, respectively) with MNE =1.0 and MNE =0.8 for L1.5^* and L2.0, respectively. These results suggest an inhomogeneous vertical distribution of aerosols, i.e., most likely the occurrence of a very large accumulation of aerosols in the ground layer if compared with the atmospheric column and with the MLH, similar to previous observations in Milan under very stable atmospheric conditions (Ferrero et al., 2011). The overestimation of scaled in situ aerosol properties compared with columnar aerosol properties observed in this study may be affected by some concurrent conditions: (a) the large role of traffic and of other ground emissions on aerosol absorption; (b) a persistent ground thermal inversion occasionally as low as a few hundred meters, according to radiosoundings at 12:00 UTC at the rural Po Valley site of San Pietro Capofiume; and (c) a bias in the ERA5 estimate of the MLH. These conditions mainly contribute to the significantly larger values observed in scaled ground absorption, particularly at 880 nm, where fossil fuel emissions provide the largest contribution. At remote sites, an opposite pattern was consistently found (e.g., Bergin et al., 2000; Slater and Dibb, 2004; Aryal et al., 2014; Chauvigné et al., 2016), where MLH-scaled surface in situ atmospheric extinction underestimates sun photometry observations of AOD, mainly because of aerosol hygroscopicity (in situ measurements are typically made at low RH) and aerosol layers above the MLH.

Figure 7 Absorption aerosol optical depth based on columnar and in situ observations in the IR (a), green (b) and UV (c) regions using both L1.5 and L2.0 AERONET retrievals. Bars indicate measurement uncertainty.

For an assessment of the role of the MLH and of atmospheric layers in the discrepancies between surface and column observations mentioned above, an analysis of the apparent aerosol optical height (ApAOH) was performed. The ApAOH can be defined as the ratio between AAOD and σ_ap (giving ApAOH units of length) to represent the atmospheric depth below which aerosols are uniformly distributed (Loría-Salazar et al., 2014). In the case of well-mixed conditions for absorbing aerosols, the ApAOH is similar to the MLH, while larger differences indicate less vertical mixing of the aerosol particles.

Figure 8 (a) Apparent aerosol optical height (ApAOH) computed based on aerosol absorption at 880 nm (870 nm) and 470 nm (440 nm) from in situ (columnar, L1.5^*) instruments. The simulated MLH depth from ERA5 is also plotted. (b) ApAOH computed based on aerosol absorption at 880 nm (870 nm) and 470 nm (440 nm) from in situ (columnar) instruments and color-coded according to the fine-mode fraction (FMF) at 440 nm. The dashed line indicates the 1:1 line.

The comparison of the ApAOH using L1.5^* and the MLH in Fig. 8a shows that ApAOH is, on average, lower than the MLH, at both the 880 (870) nm (ME $=-\mathrm{243}$ m, where ME represents the mean error) and 470 (440) nm (ME $=-\mathrm{96}$ m) wavelengths of the in situ (columnar) instruments. In the period from 22 to 26 February, the ApAOH was highly consistent with the MLH for 470 (440) nm, with ME =31 m, i.e., similar to the bias reported for ERA5 estimates of the MLH in Europe by Guo et al. (2021). Conversely, the ApAOH at 880 (870) nm was significantly lower than the MLH (ME $=-\mathrm{190}$ m), suggesting different vertical mixing between two absorbing aerosol species. The end of February 2021 featured the development of a strong anticyclone system in the Mediterranean Basin, leading to an above-average atmospheric temperature in southern and central Europe, clear sky (as shown by the high frequency of retrievals), the buildup of atmospheric pollutants and the arrival of Saharan dust-rich air masses in Italy. During the development of this high-pressure system, daily soundings collected at 00:00 and 12:00 UTC at the rural site of San Pietro Capofiume (44.65^∘ N, 11.62^∘ E; 60 km east of Modena) show the progressive vertical drop of a thermal inversion from about 2 km (on 20 February) to approximately a few hundred meters (on 26 February), leading to levels above the seasonal median for the in situ σ_ap at both 880 and 470 nm (Fig. S6). Concurrent AAOD observations were above the median at 440 nm and within the seasonal median at 870 nm, leading to different ApAOH values for these two wavelengths (Fig. 8a). Differences in the AAOD might originate from different atmospheric layering and mixing of aerosols species: the volume size distribution from the AERONET inversion shows a switch from a major modal peak in submicron diameters on 20 February (Fig. S7a) to a modal peak in the supermicron diameter range on 23 February (Fig. S7b), which lasted until 27 February at midday. Consistently, the AOD at 500 nm derived by the AERONET Spectral Deconvolution Algorithm (Sinyuk et al., 2020) had a monthly minimum during this clear-sky period, showing a switch from a fine-aerosol-controlled AOD (until 20 February) to a coarse-aerosol-controlled AOD (since 23 February). This suggests that the ApAOH at 470 (440) nm is similar to the MLH depth, most likely because of good vertical mixing of dust aerosol, as shown by the decrease in the PM${}_{\mathrm{2.5}}/$PM₁₀ ratio over the same period (Fig. S1), while the low ApAOH in the IR is probably due to the dominant contribution of (ground-level) traffic emissions to the σ_ap at this latter wavelength. This suggests that, in this case, the radiative effect due to traffic emissions was relevant mainly at the urban scale. Figure 8b shows that the ApAOH values at 880 (870) and 470 (440) nm are correlated and mainly lay on the 1:1 line during high-FMF₄₄₀ conditions and periods with an aerosol volume size distribution with a fine-mode peak (Fig. S7a). For FMF₄₄₀>0.8, a linear Pearson correlation of r=0.91 is found, with FMF₄₄₀ indicating the contribution of fine aerosol to the AOD in the blue range, where dust is a significant absorber. Under low-FMF₄₄₀ conditions, the ApAOH at 470 (440) nm increases significantly, in contrast with the ApAOH at 880 (870) nm; nonetheless, the correlation between the two remains.

3.3 Comparison of the contribution by BrC to absorption based on in situ and columnar data

The contribution of BrC to absorption in Modena according to the in situ and columnar L1.5^* data was also compared (Fig. 9). Days with a significant dust load were removed from the comparison, because the MWAA apportionment method does not include dust absorption. This is necessary because, as shown by the ApAOH, the vertical mixing of dust can be significant, affecting both columnar and in situ observations. Figure 9 compares the contribution of biomass burning to absorption in the blue region (470 and 440 nm for the columnar and in situ observations, respectively) for the two apportionment models for 17 matched data points. Calculated statistics indicate ME =0.04, MNE =1.23 and a low linear correlation (Pearson's r=0.39, Spearman's ρ=0.34). According to columnar retrievals, the biomass burning contribution ranged between 2 % and 23 %, with a median (median absolute deviation) of 3.9 % (0.8 %) and a mean (standard deviation) of 5.3 % (4.3 %); in situ observations exhibited a similar range (0 %–24 %) but suggested a higher contribution of biomass related to a lower MLH depth and a median (median absolute deviation) of 8.7 % (8.7 %). Despite the uncertainty associated with these estimates, these results highlight that urban BrC emissions have a large impact on the lower levels of the atmosphere, similar to the findings of Ferrero et al. (2011) for BC. This is consistent with the dynamics of biomass burning emissions from domestic heating, featuring a low exit velocity and negligible plume rise, particularly for natural-convection fireplaces or traditional woodstoves. This BrC absorption contribution is based on days with a negligible dust content and, thus, represents a higher-end estimate; nonetheless, it is lower than values found for polluted urban sites in eastern Asia, e.g., Beijing, Hong Kong, Seoul and Osaka, where the share of AAOD due to BrC ranged between 12 % and 14 % in the UV during non-dust days on a yearly basis (Cho et al., 2019). Even larger contributions than those found for the eastern Asian sites were reported for Europe during winter: a mean 21 % contribution of BrC to the AAOD in the blue region was reported by X. Wang et al. (2016), based on 10 years (2005–2014) of AERONET data, also excluding “dust days”. It is worth noting that, over the same decade, X. Wang et al. (2016) also reported the in situ mean share of absorption by BrC in winter at Ispra and SIRTA (Paris, France) to be ∼ 23 %. A similar study in California found that the contribution of BC and BrC to absorption by in situ and ground-based sun/sky photometers can be similar, depending on the vertical mixing of the planetary boundary layer (Chen et al., 2019). They showed that, according to both methods, the share of BrC absorption at 440 nm was approximately 30 %. Similarly, during high-pollution events in the Kathmandu Valley, Kim et al. (2021) found a good correlation between the in situ and the columnar estimates of BC and BrC. They estimated a similar share of UV absorption by BrC: 34 % and 31 % by an AE33 Aethalometer and AERONET, respectively, and a good correlation between the two techniques (R²=0.71).

Figure 9 Scatterplot of the share of total AAOD due to BrC and the share total in situ absorption due to BrC. Bars indicate measurement uncertainty.

3.4 Spatial variability in ground-based columnar retrievals

In addition to checking the vertical mixing of aerosol, the spatial mixing of ground-based columnar properties across part of the Po Valley during the study period was also investigated by comparing AERONET L1.5^* version 3 retrievals from Modena and Ispra. The observations from the two instruments were matched when collected within the same hour. The AOD data are partly scattered along the 1:1 line (Fig. 10), with a significant (at the 95 %) Pearson's (Spearman's) correlation coefficient ranging from r=0.67 (ρ=0.56) at 440 nm to r=0.84 (ρ=0.79) at 1020 nm and an orthogonal regression coefficient ranging between 0.74 at 440 nm and 1.4 at 1020 nm. The high correlation in the AOD at the two sites is mainly driven by the observations in February, when both sites experienced a drop in the FMF and an increase in the AOD because of the dust transport event. This dust event was also detected in Ispra by a ground-based lidar within the European Aerosol Research Lidar Network (EARLINET), whose total attenuated backscatter at 500 and 1064 nm showed an aerosol layer between 1.3 and 2.3 km above the ground during 23 and 24 February. The layer subsequently dropped to heights between few hundred meters above the ground and 1.8 km above the ground on 25 February (no lidar data are available for 26 and 27 February). The temporal variability in the AOD₄₄₀ and FMF₄₄₀ between the two sites also shows a similar pattern (Fig. 10). The AOD in Ispra was similar to Modena for most retrievals, with Ispra and Modena having overall median AOD₅₀₀ values of 0.27 and 0.24, respectively. The two sites differed during mid-January, when Modena had a median AOD₅₀₀=0.21, whereas Ispra had a median AOD₅₀₀=0.12, consistent with Modena being a site that is more representative of urban pollution than Ispra.

Concurrent AAOD retrievals at two sites had a different correlation pattern compared with the AOD (Fig. S8), exhibiting the largest linear correlation at 440 nm (r=0.70 and ρ=0.72, both statistically significant at the 95 % level). The correlation decreases noticeably with increasing wavelength: there is a significant correlation for 675 nm (r=0.42, ρ=0.43, 95 % significance) and 870 nm (r=0.37, ρ=0.40, 95 % significance), whereas no significant correlation occurred at 1020 nm. In contrast to the AOD, the AAOD at Ispra was larger than at Modena, with the median AAOD₄₄₀ being 0.025 and 0.011 at Ispra and Modena, respectively. The differences in the AAOD at the two sites decreased at longer wavelengths, with the median AAOD₈₇₀ being 0.003 and 0.002 at Ispra and Modena, respectively. According to the apportionment model (Fig. S5), the larger AAOD₄₄₀ at Ispra is due to the larger impact of dust at this site compared with Modena, with AAOD_440,dust being 0.022 and 0.001 at the two sites, respectively. AAOD_440,dust was the component with the largest correlation between the two sites (r=0.62, ρ=0.60, 95 % significance), followed by BC (r=0.50, ρ=0.43, 95 % significance), with similar correlation values for AAOD_870,dust and AAOD_870,BC for these two species (Fig. S5), whereas the correlation for BrC absorption at 440 nm was not significant.

Figure 10 (a–d) Comparison of hourly median AOD retrieved at Modena and Ispra at four wavelengths during the investigated period (January 2020–March 2021), color-coded according to the fine-mode fraction (FMF) at 440 nm. The dashed line indicates the 1:1 line. (e–h) Time series of the AAOD at Modena and Ispra during the investigated period (January 2020–March 2021). The bars indicate the hourly interquartile range.

4 Conclusions

In the urban area of Modena, a town representative of several urban regions in the Po Valley (a pollution hotspot for Europe), a set of light-absorbing aerosol (LAA) observations at multiple wavelengths was collected, along with meteorological and vehicle traffic data. Aerosol absorption was monitored by two in situ micro-Aethalometer instruments (MA200) and by a Cimel CE318 sun/sky photometer contributing to the AERONET network. The MA200 instruments were deployed at two locations representative of urban background and urban traffic conditions in Modena, while the Cimel sun photometer was located in an area representative of urban background conditions in Modena. In situ observations, apportioned to fossil fuel and biomass burning, were shown to be largely influenced by ground emissions. The comparison of columnar absorption and in situ absorption rescaled over the mixing layer exhibited contrasting results, demonstrated by a large difference in the infrared region (mean normalized error, MNE, of up to 2.2) but with better agreement in the blue wavelength region (MNE =0.8), confirming the impact of ground emissions on atmospheric levels of LAAs. Under the (reasonable) assumption of the generation of most of the AOD signal within the mixing layer, the heating rate due to LAAs was estimated to be 1.61 K d⁻¹. The apportionment of columnar absorption to black carbon (BC), brown carbon (BrC) and dust, along with the aerosol size distribution by AERONET inversion, highlighted the major role of long-range-transported dust in driving the correspondence between the in situ and the columnar absorption at 440 nm, indicating a deeper vertical mixing for dust, in contrast with urban ground-based emissions that are confined to lower heights. This latter result was shown specifically for BrC absorption, whose contribution to in situ absorption was larger (up to 23 %) and featured a wider variability, relative to the columnar retrieval of absorption. The spatial extent of the dust impact was evaluated by the combined analysis of concurrent columnar retrievals at Modena and Ispra (225 km northwest of Modena): the sites showed large agreement with respect to the AOD (Pearson's linear correlation coefficient of r=0.84 at 1020 nm) and AAOD at 440 nm (r=0.70), where dust has a significant absorption. Consistently, the AAOD apportioned to dust was the species with the largest correlation between the two sites, reaching r=0.62 at 440nm, supporting the occurrence of significant spatial mixing of the transported dust, along with the vertical mixing.

An improved knowledge of the role of urban emissions of LAAs is critical to control local air quality, urban heat island effects and climate forcing, and an apportionment of LAAs based on their atmospheric levels, as presented here, contributes towards this goal. This study provides important insights into the role of in situ absorption monitoring in estimating the actual absorption aloft and whether it can be used for radiative forcing estimates. Moreover, the characterization of the intra-urban variation in absorbing aerosol based on different site types contributes to the ambient exposure domain. Towards this latter outcome, a more in-depth investigation of the contribution of urban areas to atmospheric LAAs can be gained by the application of specific atmospheric dispersion tools, and this represents one of the major study outlooks. More specifically, Lagrangian particle dispersion models would provide information on atmospheric levels across the urban area at a fine spatial resolution, supporting advanced exposure studies, and, further, would give an estimate of the spatially and temporally resolved emission factors for LAAs in the urban area.

Appendix A: AAOD apportionment model

AAOD was apportioned by solving the following equation system:

$$\begin{array}{}\text{(A1)}& {\text{AAOD}}_{\mathit{\lambda }}=\sum _i{\text{AAOD}}_{\mathit{\lambda }}^i,\end{array}$$

$$\begin{array}{}\text{(A2)}& {\text{AAOD}}_{\mathit{\lambda }}^i={\text{AAOD}}_{\mathrm{ref}}{\left(\mathit{\lambda }/{\mathit{\lambda }}_{\mathrm{ref}}\right)}^{-{\text{AAE}}_{\mathit{\lambda }}^i}.\end{array}$$

Here, i is BC, dust or BrC and λ is 440, 675, or 880 nm (BrC only at 440 and 675 nm).

Table A1 Statistical metrics for the assessment of the agreement between the in situ and columnar data (L^is and $L_{\mathrm{i}}^{\mathrm{col}}$, respectively).

Appendix B: Description of symbols and acronyms used in the text

Symbol	Description
σap	Aerosol particle absorption coefficient
${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{ff}}$	σap from BC by fossil fuel combustion
${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BC},\mathrm{bb}}$	σap from BC by biomass burning
${\mathit{\sigma }}_{\mathrm{ap}}^{\mathrm{BrC}}$	σap by brown carbon
AAEc	Absorption Ångström exponent for column observations
AAEc1	AAEc across the 440–675 nm range
AAEc2	AAEc across the 675–880 nm range
AAEi	Absorption Ångström exponent for in situ observations
AAE${}_{\mathrm{ff},\mathrm{BC}}^{\mathrm{i}}$	AAEi from BC by fossil fuel combustion
AAE${}_{\mathrm{BrC}}^{\mathrm{i}}$	AAEi from BrC
AAODλ,species	Absorption AOD at wavelength λ and species (i.e., BC, BrC or dust)
ADRE	Absorptive direct radiative effect
AERONET	AErosol RObotic NETwork
AODλ	AOD at wavelength λ
ApAOH	Apparent aerosol optical height
BOA	Bottom of the atmosphere
BC	Black carbon
BL	Boundary layer
BrC	Brown carbon
DRE	Direct radiative effect
eBC	Equivalent black carbon

FMFλ	Fine-mode fraction at wavelength λ
HR	Heating rate
IQR	Interquartile range
LAA	Light-absorbing aerosol
ME	Mean error
MLH	Mixing layer height
MNE	Mean normalized error
MWAA	Multi-Wavelength Absorption Analyzer
SAEc	Scattering Ångström exponent for columnar observations
SAEc1	SAEc across the 440–675 nm range
SAEc2	SAEc across the 675–880 nm range
SAOD	Scattering AOD
SSA	Column single-scattering albedo
TOA	Top of the atmosphere
UB	Urban background
UT	Urban traffic