Estimation of PM2.5 from Satellite Observations with a Physical-Informed Method

1. Introduction

Health effects of aerosol pollution are well documented in literature [1,2], with effects ranging from allergies to chronic respiratory diseases and even carcinogenic effects as stated by WHO in 2013 [2]. For this reason, the continuous monitoring of aerosol concentrations within urban areas and in background and remote locations is of primary importance for studies focused on exposure and risk assessment practices.

Aerosol monitoring is usually performed with conventional and standardized procedures that employ daily in-situ measurements of Particulate Matter mass concentration (PM, in $\mu g/m^3$) at fixed stations. According to Martin et al. [2], in Europe there are 3 to 4 PM monitors per million of inhabitants, while in other continents the situation is even worse. As a result, most surface areas are underrepresented in terms of spatial and temporal characterization of PM.

In the last decades, the scientific community has pursued the idea of using satellite observations to increase the time and space resolutions of the sparse ground-based network [2]. It has been demonstrated that the satellite Aerosol Optical Depth (AOD) can be correlated to the ground-based PM by an empirical relation that considers meteorological variables such as the Relative Humidity (RH) and the Planetary Boundary Layer Height (PBLH). This relation can be approximated as being linear in particular cases such as in urban areas, where the aerosols are well-mixed inside the PBL and in absence of thick aerosol layers at higher altitudes [2,2]. For example, Wang et al. [2] found a correlation coefficient equal to 0.7 (out of 1) between MODIS AOD and hourly PM2.5; the correlation increased to values higher than 0.9 when lowering the time resolution to monthly averages. Koelemeijer et al. [2] considered several European cities and found a linear correlation coefficient of 0.6 between MODIS AOD and daily PM2.5, over rural and background urban stations.

More complex methodologies employ data-driven models (statistical or Machine Learning) or chemistry transport models [2]. Each methodology has its strengths and weaknesses. For example, data-driven models require huge datasets for the training phase, and it is complex to select the most significant variables to being included. On the other hand, the use of chemical transport models requires a heavy computational power, and their accuracy depends on the quality of input data such as emission sources, land use, meteorology, etc. Lastly, there are also semi-empirical methodologies, which employ measurements at selected locations to estimate the coefficients of the linear function that relates AOD and PM. This type of methodology presents low computational costs and is based on physical relationships, but the input data may be spatially and temporally sparse [9]. Nonetheless, a semi-empirical method can be implemented with some caution and tailored to the case study, depending on the available data. Ferrero et al. [10] considered factors such as the spatial representativeness of the ground-based PM with respect to the satellite footprint and the temporal difference between the satellite overpass time and the daily or hourly PM measurements, as well as the vertical characteristics of the atmosphere.

In the revised Ambient Air Quality Directive from the European Commission [11], the annual limit for $P{M}_{2.5}$ has been lowered to 10 $\mu g/m^3$, while the 24-hour average limit is 25 $\mu g/m^3$, to be exceeded less than 18 times during the year. The Po Valley, in Northern Italy, is an area known for frequent exceedances of air pollution limit values, including those for aerosol pollution. This is due to its high urbanization and industrialization density, and is exacerbated especially in the winter period by its peculiar geography bounded by the Alps and the Apennines mountains limiting atmospheric circulation, as well as characteristic weather patterns characterized by stability conditions that often trap and accumulate pollutants for long periods.

In the last decades, MODIS AOD has been used with empirical, statistical and machine learning methodologies to estimate the PM spatial variability in the area. Koelemeijer et al. [2] and Di Antonio et al. [12] managed to highlight the contrast between the high atmospheric pollution in the Po Valley compared to nearby areas. Di Nicolantonio et al. [13] was one of the first studies to implement a semi-empirical methodology, and found a correlation of 0.68 and 0.59 between MODIS AOD and daily PM2.5, reaching values of 0.7 when using monthly averages. More recent studies [10,14,15] exploited MAIAC AOD product [16], retrieved by MODIS observations with a spatial resolution of 1 km x 1 km, to better characterize the spatial variation of PM.

Among recent studies in the area, Ferrero et al. [10] and Arvani et al. [14] considered the use of a hygroscopic growth function within their methodologies and found that the results of linear correlation between AOD and PM were not significantly improved. Hygroscopicity is the property of some aerosols to grow in size according to their water uptake, changing their optical properties as well as their mass concentration. A corrective function is needed when relating the AOD, which is retrieved at ambient humidity, and the PM, which is measured in controlled environments with RH < 50%. The use of this function should improve the correlation between the two quantities. The main issue is the empirical nature of this function, whose coefficients depend on the local aerosol characteristics [17,18].

In recent years, Optical Particle Counters (OPCs) have been employed as part of air quality monitoring networks because of their portability, higher temporal resolution and limited costs with respect to conventional reference instrumentation utilized for PM measurement [19]. OPCs detect the Particle Number concentration (PN) in specific size bins and retrieve the PM mass concentration through an internal algorithm. If the OPC is operated at ambient conditions, it estimates the ambient or “wet” PM, which is supposedly better correlated with AOD than the “dry” PM from standard measurements. In fact, the hygroscopic function is not needed when relating “wet” PM with AOD at ambient conditions. OPC data has also been used as auxiliary data in previous studies such as [10], but to our knowledge, no previous study directly relates OPC-derived PM to satellite AOD.

The present study compares estimates of PM retrieved by a specific Optical Particle Counter, the LOAC (Light Optical Aerosol Counter) operated at ambient conditions, with MAIAC AOD for the year 2023 over the city of Bologna, in the Po Valley. The study also attempts to find a semi-empirical linear function to relate the "wet" PM2.5 and the satellite AOD, with the use of auxiliary data.

PN measurements from the LOAC are used as an estimate of the aerosol Particle Size Distribution (PSD), allowing the computation of the Effective Radius of the aerosol mixture. The Effective Radius is an important variable in the correlation between the optical properties and the mass concentration of aerosols and strongly affects the radiative properties derived from satellite. In addition, the LOAC infers information concerning the aerosol type (such as Salt, Mineral, and Carbonaceous aerosols) [20] which is used in the semi-empirical methodology together with Effective Radius and PBLH.

2. Materials and Methods

2.1. Study Area

In this work, ground-based observations of meteorological variables and aerosols physical characteristics were obtained during an intensive experimental field campaign conducted in 2023 in the urban area of Bologna, Italy, located in the Southern border of the Po Valley.

The Po Valley is characterized by a distinct seasonality of aerosol mass concentration. During winter, the frequent stability conditions and low wind velocity tend to favour the trapping and accumulation of pollutants [21]. Particle accumulation give rise to a concentration gradient at the PBL level, which correlates with a vertical gradient of Relative Humidity [22]. The concentration gradient is also linked to a vertical gradient of aerosol extinction. During winter, aerosol extinction tends to increase in the lowest levels of the atmosphere and concurrently to decrease at altitudes higher than 1 km; during spring and summer, the extinction is evenly distributed in the atmospheric column. As a result, during winter more than 80% of AOD refers to particles below 1 km [23].

The relative composition of PM also varies during the year according to changes in aerosol emissions. The 75% of urban PM2.5 in the Po Valley is composed of water-soluble inorganic compounds (Nitrates and Sulphates) and Carbon [24]. Local emission sources in cities comprise different types of fuel due to traffic and residential heating, which leads to different fractions of Black Carbon and organic compounds during the year [25]. In the city of Bologna, at least 40% of PM2.5 is due to traffic emissions, 24% is composed of road dust, and 8% is due to biomass burning; the levels of mineral dust, which are usually negligible, often increase in the summer during episodes of desert dust intrusion [26]. A change in the composition of the aerosol mixture affects the average optical properties retrieved from remote sensing observations, such as satellite AOD. For this reason, when comparing the aerosol mass concentration and the total extinction, it is important to obtain information about the chemical composition or the optical behaviour of the aerosol mixture.

2.2. Methodology

Based on theoretical considerations, a linear correlation between aerosol extinction and PM2.5 can be established (e.g., [2]). Formally:

$$k_{ext,aer}(z=0)=P{M}_{2.5}{\displaystyle \frac{3〈Q_{ext,aer}〉}{4{\rho }_{aer}{r}_{eff}}}$$

(1)

Equation 1 links the aerosol extinction coefficient $k_{ext,aer}$ at ground level with the corresponding mass concentration $P{M}_{2.5}$. Usually, the mean extinction coefficient and the mean extinction efficiency $〈Q_{ext,aer}〉$ are computed at visible wavelengths. All the variables are referred to ambient conditions. ${\rho }_{aer}$ is the aerosol average mass density and $r_{eff}$ is the Effective Radius of the aerosols. The Effective Radius of the PSD $n\left(r\right)$ is defined as:

$$r_{eff}={\displaystyle \frac{{\int }_0^{\infty }\pi r^{3}n\left(r\right)dr}{{\int }_0^{\infty }\pi r^{2}n\left(r\right)dr}}$$

(2)

and the average Extinction Efficiency is defined as:

$$〈Q_{ext,aer}〉={\displaystyle \frac{{\int }_0^{\infty }{Q}_{aer}\left(r\right)\pi r^{2}n\left(r\right)dr}{{\int }_0^{\infty }\pi r^{2}n\left(r\right)dr}}$$

(3)

The total AOD is computed as the integral of the aerosol extinction coefficient along the atmospheric column. In the assumption that all aerosols in the atmospheric column are well-mixed and placed below the PBL height, the vertical profile can be approximated as:

$$AOD={\int }_0^{TOA}{k}_{ext,aer}\left(z\right)dz\cong k_{ext,aer}\left(0\right)PBLH$$

(4)

so that an approximate formula relating PM2.5 to AOD (at ambient conditions) is available:

$$P{M}_{2.5,wet}\cong {\displaystyle \frac{AOD}{PBLH}}·{\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer}〉}}$$

(5)

where we underline the fact that PM is at ambient RH. The term $A={\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer}〉}}$ represents the ratio between the mass concentration and the extinction coefficient of the aerosols. Note that all variables that define the term A can be derived from OPC measurements and are therefore consistent with the measurement of $P{M}_{2.5,wet}$.

A linear regression can be performed by applying the function:

$$Y=a·X+b$$

(6)

with $Y=P{M}_{2.5,wet}$ and $X={\displaystyle \frac{AOD}{PBLH}}$. Theoretically, the slope values will vary throughout the year according to the varying aerosol mixtures and the meteorological conditions impacting on the location of interest. It has been verified that, in general, the higher the PBLH, for example during summer, the lower the linear correlation (e.g. [27]). When using Eq. 6 to derive the empirical coefficient, the slope is considered to be constant in the time window and location considered by the dataset. Lastly, b is the offset between $P{M}_{2.5,wet}$ and $AOD$ and, in theory, it will not vary seasonally.

The present study also analyses the linear regression between PM2.5 and the AOD corrected by all available auxiliary data: Effective Radius, average mass density and average optical properties, thus having an estimate of the term A. The linear correlation is performed with Eq. 6 and with $Y=P{M}_{2.5,wet}$ and alternatively with $X={\displaystyle \frac{AOD}{PBLH}}$ and $AX=A{\displaystyle \frac{AOD}{PBLH}}={\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer}〉}}{\displaystyle \frac{AOD}{PBLH}}$ and the results are compared.

PM2.5 absolute uncertainty is fixed as the standard deviation and AOD relative uncertainty is 0.05 (see Sec. 2.4). Uncertainties for $X_1={\displaystyle \frac{AOD}{PBLH}}$ and $X_2=A{X}_1$ are computed with the relative uncertainty rule:

$${\displaystyle \frac{d{X}_1}{X_1}}=\sqrt{0.{05}^2+{\left({\displaystyle \frac{dPBLH}{PBLH}}\right)}^2}{\displaystyle \frac{d{X}_2}{X_2}}=\sqrt{{\left({\displaystyle \frac{d{X}_1}{X_1}}\right)}^2+{\left({\displaystyle \frac{d{r}_{eff}}{r_{eff}}}\right)}^2}$$

(7)

2.3. OPC Data

OPC measurements were obtained by a LOAC (Light Optical Aerosol Counter) [20] mounted on the rooftop of the Department of Physics and Astronomy (DIFA) of the University of Bologna (lat: 44.499, lon: 11.354), at about 22 m from the ground (from now on named only DIFA). The surrounding area is a vegetated street inside the historical centre of Bologna (Figure 1). The LOAC dataset covers the months from February to September 2023. The LOAC retrieves the Particle Number (PN) concentration ($\#/c{m}^3$) in 19 size bins for particle diameters between 0.2 and 30 $\mu m$ (Table 1) by the use of laser beams at 650 nm wavelength. By using two scattering angles (12° and 60°) it also determines the degree of absorption by the detected particle, from transparent to highly absorbing. This estimate is used to classify the particles, divided in five size classes (Table 2), in four types: “droplet”, “salt”, “mineral”, “carbon”. From the PN measurements and the particle type, the PM1, PM2.5, and PM10 ($\mu g/m^3$) are calculated with an internal algorithm which uses a suitable value of particle mass density ($g/c{m}^3$) according to particle type. In general, all OPC instruments are highly affected by high Relative Humidity, cloudy conditions and precipitation events (eg. [28]). As such, recent studies recommend to perform measurements with such sensors under low time resolutions (hourly and daily averages) and clean sky conditions [29].

The 1-minute PN and PM2.5 data is screened with a Hampel filter on a 10-minute window [30] to remove outliers and then averaged hourly. The Effective Radius $r_{eff}$ is also computed by discretizing Eq. 2, which becomes:

$$r_{eff}\cong {\displaystyle \frac{{\sum }_{i}P{N}_i{R}_i^3}{{\sum }_{i}P{N}_i{R}_i^2}}$$

(8)

where $P{N}_i$ is the Particle Number concentration in the size bin i and $R_i$ is half the nominal diameter of size bin i. This is because the Particle Number concentration in size bin i is the total value of $n\left(r\right)$ within the size bin: $P{N}_i\cong {\int }_{D_i/2}^{D_{i+1}/2}n\left(r\right)dr$. The $r_{eff}$ is computed with the original filtered 1-minute PN data and then averaged hourly.

The particle type dataset from LOAC is used together with PN data to compute the average aerosol Extinction Efficiency and mass density. This is done by computing

• the Number Ratio (NR) of each particle type in each size bin,
• the average single particle Extinction Efficiency of each particle type in each size bin

and by choosing the proper mass density for each particle type, assuming it is the same in all size bins. In Renard et al. 2016 [20] it is advised to use particle type information with a statistical approach. For each hour, we consider the number of occurrences of each particle type in each size bin, divide it by the total number of observations in each size bin in that hour (maximum of 60), and convert the fraction into decimal percentage. In this way, we define the Number Ratio $F_i^x$ in size bin i for the particle type x.

In Renard et al. 2010 [31], the absorption of different particle types is defined in relation to different values of the imaginary part k of the Refractive Index (RI). The “salt” and “mineral” particle type match simulations of optical properties for glass beads with $k=0.005$, while the “carbon” particle type match simulations with $k=0.43$, considered at the laser wavelength of 650 nm. We consider the optical properties of aerosols as defined in the database OPAC [32]. The urban aerosols are defined as a mixture of Soluble and Insoluble aerosols, characterised by low absorption at 650 nm ($k=0.004$ and $k=0.008$) and of Black Carbon (or Soot) which has a high absorption at 650 nm ($k=0.4$). In this study, we assume that each LOAC particle type can be related to each OPAC aerosol type:

• “salt” LOAC type related to “Soluble” OPAC type (here considered at RH = 50%),
• “mineral” LOAC type related to “Insoluble” OPAC type,
• “carbon” LOAC type related to “Soot” OPAC type,

and use the related single particle Extinction Efficiency $Q_{aer}\left(r\right)$ values.

$Q_{aer}\left(r\right)$ varies with particle radius r in different ways for each particle. We compute $Q_i^x$ as the average value of $Q_{aer}$ in the size bin i for type x. Lastly, the average Extinction Efficiency for each LOAC observation is discretized from Eq. 3 as:

$$〈Q_{ext,aer}〉\cong {\displaystyle \frac{{\sum }_i{Q}_i^1{F}_i^{1}P{N}_i+Q_i^2{F}_i^{2}P{N}_i+Q_i^3{F}_i^{3}P{N}_i}{{\sum }_{i}P{N}_i{R}_i^2}}$$

(9)

while the average mass density ${\rho }_{aer}$ is computed as:

$${\rho }_{aer}\cong {\displaystyle \frac{{\sum }_{i}P{N}_i(m_i^1{F}_i^1+m_i^2{F}_i^2+m_i^3{F}_i^3)}{{\sum }_{i}P{N}_i}}$$

(10)

where $m_i^x$ are the nominal values of mass density for each LOAC particle type x (same value for all size bins). We choose to use the same mass density values used in the computation of the PM instead of the mass density values reported in OPAC for the chosen aerosol types.

2.4. AOD

MAIAC AOD is retrieved at 550 nm from MODIS observations at high spatial resolution (1 km x 1 km) on a fixed grid. It has a better performance over bright surfaces compared to other algorithms and has been previously applied over urban areas [33]. In cases of anthropogenic aerosols, AOD below 0.12 can be overestimated. AOD uncertainty indicated in the product files is related to surface brightness. A value below 0.05 allows the algorithm to perform a standard retrieval [33,34]. A 5% uncertainty value for the single retrieval values can be used. More information on the retrieval can be found at [16]. Intense pollution cases in urban areas within the Po Valley, especially in very high polluted cities such as Milan and Turin, can be sometimes interpreted as clouds and filtered out [14].

Co-location of AOD and ground-based data at DIFA site is performed on the fixed grid used in MAIAC retrievals. The cell containing the DIFA location is selected (Figure 1). This grid cell refers to the innermost part of Bologna historical centre, characterised by medium-bright surfaces, traffic and low industrial emissions. The AOD data is acquired for the whole period covered by the LOAC dataset (February-September 2023). The total number of satellite overpasses in the study period is 985, but only 26.8% is characterized by a successful AOD retrieval. The available retrievals are few during April and May (15-16%) and more frequent in August (41%).

Considering the dataset of co-located AOD retrievals, AOD uncertainty due to surface brightness is lower than 0.045 with a median of 0.016, and AOD values are lower than 0.45 with a median of 0.14. This means that the retrievals used in the following analysis are computed with the standard algorithm, and likely not affected by underestimation.

2.5. Auxiliary Variables

Planetary Boundary Layer Height (PBLH) is retrieved by a ceilometer Vaisala CL31, mounted on the rooftop of DIFA alongside the LOAC. The ceilometer continuously measures the vertical backscatter signal from the ground level to a maximum of 7.5 km and provides estimates of PBLH with a temporal resolution of 16 s and a vertical resolution of 10 m. It also provides an estimate of Cloud Height if clouds are detected. The internal algorithm is based on the difference in backscatter signal between the mixing layer, where the signal is higher due to aerosols, and the free atmosphere, where the backscatter signal abruptly drops to very small values in absence of scattering layers [35]. An incorrect PBLH detection can be related to the presence of residual layers or to PBL changes during the early evening. In these cases, the ceilometer detects more than one PBLH value, with the highest one usually corresponding to residual layers or incorrectly detected clouds [36].

The original dataset is averaged on a 10 minutes time frame around the satellite overpass and the PBLH is selected from the retrieved Boundary Layer heights. Three possible PBLH values are retrieved by the Ceilometer algorithm. The lowest value usually refers to the strong gradient between the ground-level atmospheric pollution and the aerosols mixed above, or to a consistent residual layer (especially in winter), with values below 600 m. The middle layer, when present, is the most probable PBLH value if below 1.5 km, otherwise it may refer to clouds during precipitation or cloudy days. The third layer usually refers to clouds or extinction gradients due to complex vertical profiles. In this work, we choose the first PBLH value if higher than 600 m and the second one if below 2 km. This is consistent with the average PBLH values in the area (as mentioned for example in [23]).

Auxiliary data include quality-proofed data from the local environmental agency (ARPAE) [37]: 30-minute Relative Humidity and hourly Precipitation Rate are measured from the closest meteorological station, Bologna Idrografico (lat: 44.499, lon: 11.346), located 620 m from DIFA; daily "dry" PM2.5 from ARPAE Porta San Felice station (lat: 44.499, lon: 11.328), located 2120 m from DIFA. In Figure 1 the stations are shown alongside DIFA and the AOD cell from satellite.

Lastly, we use data from AERONET [38], a global network of ground-based sunphotometers that measure the AOD and related variables [39]. AERONET nearest stations are in the urban area of Modena (lat: 44.629, lon: 10.948) and in the rural area of San Pietro Capofiume (lat: 44.650, lon: 11.620).

3. Results

3.1. Analysis of the Particle Concentration Variability

In the following, we analyze the data obtained by the LOAC. The measurement period considered in this work (February-September 2023) was characterized by a particular variability in the meteorological conditions impacting on the study area.

The month of May 2023 was affected by three exceptional rainfall events (May 2, 16, and 17), with cumulative daily rainfall values higher than 100 mm recorded at the Bologna Idrografico station. Other minor rainfall events occurred throughout the measurement period. Figure 2 shows daily averages of “wet” PM2.5 from LOAC and RH from Bologna Idrografico for the whole period. Precipitation events with rainfall higher than 20 mm are highlighted in grey. As noted in [29], LOAC measurements are affected by the hygroscopic growth of particles under precipitation events and high RH conditions. This is apparent from the daily values of PM2.5 recorded during the study period, which are all below 35 $\mu g/m^3$ except for observation records during precipitation events, in which PM2.5 reaches values above 50 $\mu g/m^3$.

Desert dust intrusions were frequent during the months of June and July. The monthly PSD averages in Figure 3 shows a corresponding increase in particles with radius greater than $1\mu m$, compared to particles with smaller radii, during these months. PSD is computed by dividing the PN by the width of each size bin, as explained in Section 2. The Effective Radius, computed with Eq. 8, shows values larger than $2\mu m$ during precipitation events in May 2 and 16, and higher than $1.5\mu m$ between June $19-22$, and July $9-11$ and $15-18$. Data from AERONET nearest stations of Modena and San Pietro Capofiume highlight that, during that time, the Angstrom Exponent $\alpha$ was below 0.5 in June and below 0.6 in July. Low values of $\alpha$ are associated with desert dust transport which is frequent in Italy during summer, with dust plumes that also reach the Po Valley and Bologna [40]. Except for precipitation and dust transport events, the Effective Radius is less than $0.6\mu m$ during the measurement period.

Figure 5 shows the monthly particle type occurrences, described by the Number Ratios as defined in Section 2. Only the first three size ranges are shown, consisting of particles with diameters between $0.2$ and $5\mu m$ (as specified in Table 2). Particles in the fourth size range ($5-12.5\mu m$) consist of water droplets or unidentified types. Particle types in the two largest size ranges ($12.5-40\mu m$) are mostly not identified, probably due to measured PN concentrations lower than $0.01\#/c{m}^3$. The “salt” and “mineral” types are present in the first two size ranges in varying proportions during the year. In the third size range ($0.9-5\mu m$), “mineral” type is the most common, but also the “carbon” type is consistently present.

When calculating $〈Q_{ext,aer}〉$ and ${\rho }_{aer}$, water droplets are not considered, while the values for larger size bins are set to constant values due to the high frequency of unidentified particle types.

Figure 4. Daily Effective Radius from LOAC (blue line). Grey areas refer to days with total precipitation greater than 20 mm. Red areas refer to dust transport events.

Figure 4. Daily Effective Radius from LOAC (blue line). Grey areas refer to days with total precipitation greater than 20 mm. Red areas refer to dust transport events.

Figure 5. Monthly Number Ratio of the detected particle types in the first 3 size classes. Blue = droplet, Red = salt, Yellow = mineral, Purple = carbon. (a) Size class $0.2-0.5\mu m$. (b) Size class $0.5-0.9\mu m$. (c) Size class $0.9-5\mu m$.

Figure 5. Monthly Number Ratio of the detected particle types in the first 3 size classes. Blue = droplet, Red = salt, Yellow = mineral, Purple = carbon. (a) Size class $0.2-0.5\mu m$. (b) Size class $0.5-0.9\mu m$. (c) Size class $0.9-5\mu m$.

3.2. Methodology Application

The present section explains the assumptions and the results of the linear regressions. The goodness of the correlation is evaluated using the Pearson’s correlation coefficient R (95% confidence) and the linear regression performance is measured by the coefficient of determination $R^2$ and Root-Mean Square Error (RMSE). The linear regression function and the correlated variables are defined in Eq. 6. We also normalize the coefficient A considering the different units of measurement between PM2.5 ($\mu g/m^3$), ${\displaystyle \frac{AOD}{PBLH}}$ ($m^{-3}$), $r_{eff}$ ($\mu m$) and ${\rho }_{aer}$ ($g/c{m}^3$); by definition, $〈Q_{ext,aer}〉$ is unitless.

The uncertainty associated with the hourly averages of $P{M}_{2.5,wet}$ and $r_{eff}$ is the standard deviation. During particular events such as precipitation and dust transport, the uncertainties are higher because of the higher variability in the PN measurements. For periods when the AOD data is available, the relative uncertainty of $P{M}_{2.5,wet}$ is below 40% with an average of 36% and the relative uncertainty of $r_{eff}$ is lower than 30% with an average of 22%. Only this data is considered in the methodology because the AOD is needed to perform the linear regressions.

On the other hand, AOD has a nominal relative uncertainty of 5% and the composed uncertainty of $X={\displaystyle \frac{AOD}{PBLH}}$ is between 5% and 10% in the whole AOD-PBLH dataset. The uncertainty for $AX$ is between 5% and 50%, which can be higher than the uncertainty on Y. This is expected as it is a composed uncertainty computed from X and $r_{eff}$.

First, we perform the linear regression with $Y=P{M}_{2.5,wet}$ and $X={\displaystyle \frac{AOD}{PBLH}}$, without adding auxiliary information. We consider PM2.5 data whose standard deviation is less than $3\mu g/m^3$ and whose relative uncertainty is less than 20%. We also consider AOD data between 0.01 and 0.35 to avoid errors in the retrieval algorithm. PBLH values under 600 m are considered faulty retrievals and set to 800 m. By using these rules, May and September data are filtered out, as well as several data points during June and July. The final dataset therefore consists of a total of 56 observations available.

The Pearson’s correlation coefficient R is 0.5587 (out of 1), indicating a weak linear correlation between the two variables. The linear regression function has $RMSE=4.28\mu g/m^3$ and $R^2=0.312$, which is not enough for a good linear fit. The linear function is $Y_1=7.26+0.282·X_1$. The linear function is shown in Figure 6a with the selected data points. As expected, the relationship between the variables does not show any clear dependence on RH, suggersting that there is no need for a hygroscopic correction factor.

For comparison, we perform the linear regression using daily "dry" PM2.5 from reference ARPAE station Porta San Felice in the same temporal window. As PM2.5 data is obtained from reference instruments, we assume that the uncertainty is always below 20%. We assume that the PBLH does not change significantly between the station and DIFA, therefore using PBLH obtained from the ceilometer. The selected AOD refers to the cell whose centre is closest to the ground-based station. To obtain daily averages of $X={\displaystyle \frac{AOD}{PBLH}}$, we first compute the fraction and then average daily the available data. In this case, the regression is performed on 102 data points between February and September 2023. The value of R is 0.1558, indicating an absence of linear correlation. As it is apparent from Figure 6b, "dry" PM2.5 and the corrected AOD seem uncorrelated, and there is no correlation with the respective RH value (daily averaged). The linear function is $Y_1^{\prime }=10+0.15X_1^{\prime }$ with $R^2=0.0243$ and $RMSE=6.17\mu g/m^3$.

Now we consider LOAC PM2.5 and the auxiliary information to compute $AX$ and its uncertainty. We recall that $A={\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer}〉}}$ is the theoretical coefficient between "wet" PM2.5 and ${\displaystyle \frac{AOD}{PBLH}}$. We filter data whose uncertainty is higher than $2\mu g/m^3$. The filtered subset comprises 49 data points between February and July. The Pearson’s correlation coefficient is 0.7610, indicating a much higher linear agreement between the variables. The linear regression function has $RMSE=3.11\mu g/m^3$ and $R^2=0.579$, which indicates a slightly better linear fit compared to results from previous studies such as [13,41]. The linear function, presented in Figure 7a, is $Y=6.27+0.812·AX$.

The offsets b in both linear functions are compatible with each other within confidence bounds (see also Table 3). This means that this value depends on other factors, such as the difference between the point measurements (LOAC and ceilometer data) and the satellite AOD retrieval (with a footprint of $1k{m}^2$). Theoretically, the slope of the second linear function should be close to 1, as A represents the theoretical slope that relates PM2.5 to AOD and should account for all the variability. The value of 0.812 can be interpreted as the effects of the assumptions made for the methodology application.

It should be noted that the smallest particle diameter detected by the LOAC is $0.2\mu m$ (Table 1). This excludes the detection of particles with radii less than $0.1\mu m$, that are almost negligible in the computation of mass concentration, but can have a significant contribution in terms of particle number concentration and in the extinction of radiation. Moreover, the Effective Radius should be calculated considering particles with radii starting from $0.05\mu m$, for example, as in the definition of aerosol types in OPAC [32]. The same is valid for the computation of $〈Q_{ext,aer}〉$ and ${\rho }_{aer}$.

To improve the calculation of these parameters, we add a smaller size bin for diameters $0.1-0.2\mu m$, whose Number Ratio for each particle type is assumed to be the same as the first size range ($0.2-0.5\mu m$) and whose Particle Number concentration is extrapolated by considering PN values in size bins 1 to 5. Taking into account only the 49 data points used for the second linear regression, $r_{eff,new}$ values fall in the $50-90\%$ range of the original $r_{eff}$, $〈Q_{ext,aer,new}〉$ values are $50-75\%$ of the original $〈Q_{ext,aer}〉$, and ${\rho }_{aer,new}$ values are nearly equal to the original ones. We repeat the linear regression with the new values of $r_{eff}$, computing $A_1={\displaystyle \frac{4{\rho }_{aer}{r}_{eff,new}}{3〈Q_{ext,aer}〉}}$. In this case, the slope of the linear function reaches a value of $a=1.042$ with $R^2=0.538$, $RMSE=3.03\mu g/m^3$ and $R=0.7337$ (shown in Figure 7b). To complete the analysis, we also compute $A_2={\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer,new}〉}}$. In this case the linear regression has a slope $a=0.431$ with $R^2=0.49$, $RMSE=3.18\mu g/m^3$ and $R=0.6997$. Finally, when using the new values for both variables, we compute $A_3={\displaystyle \frac{4{\rho }_{aer}{r}_{eff,new}}{3〈Q_{ext,aer,new}〉}}$ and in this case the slope is $a=0.587$, with $R^2=0.518$, $RMSE=3.09\mu g/m^3$ and $R=0.7196$. Table 3 summarizes the linear correlation results obtained in each case considered so far. The fit results improve when using the new computed values of $r_{eff}$ while they do not improve when considering also the new computed values of $〈Q_{ext,aer}〉$.

4. Discussion

We show that the use of ’wet’ PM2.5 mass concentration observations is correlated to the collocated satellite AOD, thus avoiding the application of hygroscopic corrections to the ground-based data. The Pearson’s correlation coefficient $R=0.56$ shows that the "wet" hourly PM2.5 is linearly correlated to AOD, even without using auxiliary information with the exception of the PBL values which is always considerd. For comparison, "dry" PM2.5 from reference ground-based measurements only yields a coefficient $R=0.16$, indicating the absence of a linear correlation with satellite AOD.

The "wet" PM2.5 is representative of the mass concentration of the particulate matter in humidity conditions, consistent with the AOD satellite measurement. The analysed data show that higher values of both PM and AOD are generally associated to higher values of RH, contrary to "dry" PM2.5 which is measured at fixed RH values lower than 50%. This allows us to use data from all year instead of selecting data by season, as the main seasonal differences depend on RH, as well as other meteorological variables such as Atmospheric Temperature and Wind Speed. This is an element of novelty with respect to previous works on the same topic, while reaching results comparable to previous studies [2,10,13,14].

Moreover, a Pearson’s correlation coefficient of $R=0.76$ is obtained when AOD is AOD is multiplied by parameters obtained by the LOAC. An increase in R points at a better linear correlation between the variables. The coefficient of determination of $R^2=0.58$ agrees with previous studies that used a linear fit for the correlation [13,41], with the addition that in our case the linear function is valid for the entire temporal window investigated (February-September 2023), without the need of using different coefficients depending on the specific season or the month.

The slope of the linear function between "wet" PM2.5 and adjusted AOD is 0.82, which is lower than the perfect fit value of 1. This means that the applied corrections do not fully relate the two variables linearly. This is an expected conclusion due to the high relative uncertainty of the computed values of the theoretical slope $A={\displaystyle \frac{4{\rho }_{aer}{r}_{eff}}{3〈Q_{ext,aer}〉}}$, which in turn is due to the high relative uncertainty of $r_{eff}$ with respect to PM2.5, and to the absence of a quantitative estimate for the uncertainty of the particle type inferred from the LOAC measurements. Moreover, the linear approximation may not hold for all possible cases.

By repeating the linear regression with an alternative calculation of $r_{eff}$, the slope of the linear function reaches a value of 1.04 while maintaining a comparable goodness of fit ($R^2=0.54$), relative error on coefficients ($24-27\%$), and Pearson’s coefficient ($R=0.73$). The alternative $r_{eff}$ is calculated by extrapolating Particle Number concentration for smaller particles, which are not measured by the LOAC, but are not negligible in the computation of aerosol optical properties. The new $r_{eff}$ values are lower than the original ones while maintaining comparable uncertainties. This suggests that using PN measurements from LOAC may overestimate the value of $r_{eff}$ in urban environments.

Our methodology takes into account the theoretical relationship between PM mass concentration and total aerosol extinction, and is able to correlate them with a slope close to 1. Similarly to previous studies, an offset is needed in the linear correlation to account for factors other than aerosol properties and meteorology (e.g., local traffic).

An appropriate correction is also needed to compare the "wet" PM2.5 retrieved from satellite data with standard measurements of “dry” PM2.5. It should consist of an empirical function whose coefficients are related and fine-tuned in accordance with the aerosols features of the considered location. We argue that the conversion from "wet" to "dry" PM2.5 is secondary to the estimation of PM from satellite data and therefore it is not furtherly studied in the present work.