Assessment of the impact of NO₂ contribution on aerosol-optical-depth measurements at several sites worldwide

This work aims at investigating the effect of NO₂ absorption on aerosol-optical-depth (AOD) measurements and Ångström exponent (AE) retrievals of sun photometers by the synergistic use of accurate NO₂ characterization for optical-depth estimation from co-located ground-based measurements. The analysis was performed for ∼ 7 years (2017–2023) at several sites worldwide for the AOD measurements and AE retrievals by Aerosol Robotic Network (AERONET) sun photometers which use OMI (Ozone Monitoring Instrument) climatology for NO₂ representation. The differences in AOD and AE retrievals by NO₂ absorption are accounted for using high-frequency columnar NO₂ measurements by a co-located Pandora spectroradiometer belonging to the Pandonia Global Network (PGN). NO₂ absorption affects the AOD measurements in UV-Vis (visible) range, and we found that the AOD bias is the most affected at 380 nm by NO₂ differences, followed by 440, 340, and 500 nm, respectively. AERONET AOD was found to be overestimated in half of the cases, while also underestimated in other cases as an impact of the NO₂ difference from “real” (PGN NO₂) values. Overestimations or underestimations are relatively low. About one-third of these stations showed a mean difference in NO₂ and AOD (at 380 and 440 nm) above 0.5 × 10⁻⁴ mol m⁻² and 0.002, respectively, which can be considered a systematic contribution to the uncertainties in the AOD measurements that are reported to be of the order of 0.01. However, under extreme NO₂ loading scenarios (i.e. 10 % highest differences) at highly urbanized/industrialized locations, even higher AOD differences were observed that were at the limit of or higher than the reported 0.01 uncertainty in the AOD measurement. PGN NO₂-based sensitivity analysis of AOD difference suggested that for PGN NO₂ varying between 2 × 10⁻⁴ and 8 × 10⁻⁴ mol m⁻², the median AOD differences were found to rise above 0.01 (even above 0.02) with the increase in NO₂ threshold (i.e. the lower limit from 2 × 10⁻⁴ to 8 × 10⁻⁴ mol m⁻²). The AOD-derivative product, AE, was also affected by the NO₂ correction (discrepancies between the AERONET OMI climatological representation of NO₂ values and the real PGN NO₂ measurements) on the spectral AOD. Normalized frequency distribution of AE (at 440–870 and 340–440 nm wavelength pair) was found to be narrower for a broader AOD distribution for some stations, and vice versa for other stations, and a higher relative error at the shorter wavelength (among the wavelength pairs used for AE estimation) led to a shift in the peak of the AE difference distribution towards a higher positive value, while a higher relative error at a lower wavelength shifted the AE difference distribution to a negative value for the AOD overestimation case, and vice versa for the AOD underestimation case. For rural locations, the mean NO₂ differences were found to be mostly below 0.50 × 10⁻⁴ mol m⁻², with the corresponding AOD differences being below 0.002, and in extreme NO₂ loading scenarios, it went above this value and reached above 1.00 × 10⁻⁴ mol m⁻² for some stations, leading to higher AOD differences but below 0.005. Finally, AOD and AE trends were calculated based on the original AERONET AOD (based on AERONET OMI climatological NO₂), and its comparison with the mean differences in the AERONET and PGN NO₂-corrected AOD was indicative of how NO₂ correction could potentially affect realistic AOD trends.

1 Introduction

Earth's radiation budget and climate is impacted by both the direct and indirect effects of atmospheric aerosols (IPCC, 2021). The direct effect of aerosols is associated with the absorption and scattering of solar radiation (Hobbs, 1993), while the indirect effect involves the interaction of aerosols with clouds by acting as cloud condensation nuclei and potentially altering cloud properties, precipitation, surface fluxes, and the energy budget of the atmosphere (Rosenfeld et al., 2014; Herbert and Stier, 2023). Apart from the impact on climate and radiative forcing, aerosols also have adverse effects on human health, leading to respiratory, cardiovascular, and neurological diseases, hypertension, diabetes, and even cancer (Lelieveld et al., 2015; Molina et al., 2020). Aerosol optical depth (AOD) is the most widely used parameter for the estimation of columnar atmospheric aerosol concentrations at different spectral bandwidths.

Sun photometers are passive remote sensing instruments that are used for measuring AOD, which is calculated using the Lambert–Beer law by taking into account the contribution from Rayleigh scattering by atmospheric molecules and absorption by atmospheric constituents other than aerosols like ozone, nitrogen dioxide, and water vapour. The global aerosol networks such as AERONET (Aerosol Robotic Network; https://aeronet.gsfc.nasa.gov, last access: 6 September 2024), SKYNET (https://www.skynet-isdc.org/aboutSKYNET.php, last access: 6 September 2024; Nakajima et al., 2020), and GAW-PFR (Global Atmospheric Watch – Precision Filter Radiometer Network; Kazadzis et al., 2018) use specific methodology to account for the optical-depth contributions from these atmospheric constituents in order to retrieve AOD.

AERONET performs optical-depth corrections for Rayleigh scattering at all wavelengths: ozone for the spectral range 340–675 nm, NO₂ for the spectral range 340–500 nm, water vapour for 1020–1640 nm, and carbon dioxide and methane for 1640 nm. The uncertainty in AOD measurements from the AERONET algorithm is estimated to be ∼ 0.01 in the visible range that reaches up to ∼ 0.02 in the UV region (Eck et al., 1999; Giles et al., 2019). Other factors contributing to the AOD uncertainty in different spectral bands include the optical-depth estimation from trace-gas (ozone and NO₂) absorption, which is sensitive to the estimation of the gas concentrations. Specifically, NO₂ absorption is predominant in lower wavelengths (340–500 nm), and hence, the NO₂ correction is of significant importance at these wavelengths. This enhances the need to investigate the impact of NO₂-absorption-based optical-depth on AOD measurements and the possibility of improvements in the retrieval algorithm with a more accurate NO₂ optical-depth estimation using ground-based NO₂ measurements.

Emissions of nitrogen oxides on a global scale from natural sources are more significant than that generated from anthropogenic activities (Seinfeld and Pandis, 2016). The natural sources of NO_x emissions include wildfires, lightning, oxidation of biogenic ammonia, and microbial processes in soils. The NO₂ levels due to NO_x emissions from natural sources are referred to as background and are smaller in magnitude in comparison to the anthropogenic NO_x emissions (Koukouli et al., 2022). The NO_x budget is dominated by fossil fuel combustion, biomass burning emissions, and anthropogenic activities.

Due to inhomogeneous local emission patterns and photochemical destruction in heavily polluted regions, the NO₂ has high spatiotemporal variations and a shorter lifetime, with regional confinement near its source (Richter et al., 2005; Boersma et al., 2008; Tzortziou et al., 2014, 2015; Drosoglou et al., 2017; Fan et al., 2021). The high spatiotemporal variation in the tropospheric NO₂ can produce significant bias in the AOD measurements (Arola and Koskela, 2004; Boersma et al., 2004). Therefore, the regions with high tropospheric NO₂ emissions will have a higher likelihood for deviation from the climatological mean values (Giles et al., 2019). Furthermore, there can also be significant diurnal variation in NO₂ (Boersma et al., 2008). Hence, the climatological mean NO₂ values might not be able to represent the actual NO₂ loading and spatial distribution in the atmosphere. This in turn tends to produce potential errors in the calculation of AOD in the spectral regions with significant NO₂ absorption. However, synergistic assistance from the models, satellite observations, or co-located surface-based measuring instruments capable of providing temporal columnar products of NO₂ can help in the reduction in the associated uncertainty, and hence, the accuracy of the total column NO₂-optical-depth estimation can increase (Herman et al., 2009; Tzortziou et al., 2012). To this end, the Pandonia Global Network (PGN) (https://www.pandonia-global-network.org, last access: 6 September 2024), which is a global network of Pandora spectroradiometers that are used for trace-gas measurements and provide the NO₂ concentration, can be useful. These instruments can be used to provide a good estimation of the NO₂ concentration in the atmosphere that can help reduce the uncertainty in AOD measurements.

Here we try to follow up on a previous work by Drosoglou et al. (2023b), who analysed the impact of NO₂ absorption using PGN-spectroradiometer-based high-frequency columnar NO₂ on AOD, Ångström exponent (AE), and single-scattering albedo (SSA) retrievals from AERONET and SKYNET for the urban area of Rome (Italy) for a time period from 2017–2022. The NO₂-based AOD correction showed a systematic overestimation of AOD and AE, with mean AOD bias of ∼ 0.003 and ∼ 0.002 at 380 and 440 nm, respectively, for AERONET, a much higher (∼ 0.007) bias for SKYNET, and an average AE bias of ∼ 0.02 and ∼ 0.05 for AERONET and SKYNET, respectively. However, for high-columnar NO₂ concentrations (> 0.7 DU, Dobson unit), the average AOD bias ranged between 0.009 and 0.012 for AERONET and ∼ 0.018 for SKYNET. As this study was limited to only one location, a worldwide analysis is needed to better analyse such an NO₂-correction-based bias in AOD measurements.

The work presented in this article deals with updating the work of Drosoglou et al. (2023b) that was based in only one station. This is a first attempt to analyse a worldwide scenario in which AERONET and PGN instruments are co-located. So a more specific investigation is performed on a worldwide scale for evaluating the effect of low- to high-NO₂ loads on the AOD measurements by ground-based remote sensing in several sites across the globe in order to understand the wider impact of uncertainties introduced in the aerosol property retrievals by NO₂ absorption. In particular, we analyse a long-term dataset (∼ 7 years) collected from 33 sites distributed worldwide and where co-located measurements of both NO₂ from Pandora spectroradiometers as part of PGN and AOD from AERONET sun photometers are available. Following the Introduction, Sect. 2 deals with the observational data, methodology for the co-located stations, and the retrieval of the aerosol parameters used for the analysis and trend analysis, followed by Sect. 3, which presents the results and discussions. Finally, Sect. 4 summarizes the findings of this study.

2 Data and methodology

2.1 Data

2.1.1 Columnar aerosol property measurements (AOD and AE)

AERONET provides the datasets of aerosol optical, microphysical, and radiative properties through ground-based passive remote sensing using Cimel sun photometers (https://www.cimel.fr/solutions/ce318-t/, last access: 6 September 2024). It has a centralized data processing and distribution system providing the instrument calibration standardization and data acquisition. AERONET direct-sun algorithm data products obtained from the version 3 processing algorithm (Giles et al., 2019) are employed in this work, including level 1.5 AOD measurements at 340, 380, 440, 500, 675, and 870 nm and AE retrievals at 440–870 and 340–440 nm. Level 1.5 data products are cloud-screened and quality-assured. AERONET data used in this work cover a time period between 2017 and 2023, during which synchronous data from the co-located PGN Pandora instrument are also available. For the trend analysis in Sect. 2.2.3, AERONET AOD data between 2013–2023 are considered. The standard AERONET AOD calculations are based on the NO₂-optical-depth estimation from the Ozone Monitoring Instrument (OMI/Aura) level 3 climatological (hereafter referred to as OMIc) total NO₂ values at a spatial resolution of 0.25° × 0.25° and for the time period between 2004 and 2013.

2.1.2 Vertical column NO₂ measurements

The total NO₂ column product used in this study is obtained from Pandora spectroradiometers which are part of PGN. Pandora spectroradiometers perform direct solar irradiance and scattered radiance measurements with a high temporal resolution in the spectral range of 280–530 nm for the retrieval of tropospheric and total column densities, near-surface concentrations, and vertical profiles of atmospheric trace gases (e.g. NO₂, O₃, and HCHO) (e.g. Herman et al., 2009; Tzortziou et al., 2012, 2015). The total column NO₂ densities are retrieved from the direct-sun measurements with ∼ 0.6 nm resolution in the spectral range of 280–530 nm, using Blick software by Cede (2024). Pandora NO₂ vertical column density (VCD) used in this analysis is obtained from level 2 datasets that provide column amounts, concentrations, profiles, a direct-sun retrieval code “nvs3”, and Blick processor version 1.8. From this dataset, total column NO₂ VCD with high- (0, 10) and medium-quality (1, 11) flags are considered.

2.1.3 Satellite observations

Daily tropospheric NO₂ columns are retrieved from OMI/Aura level 3, version 1.1, global data products gridded as 0.25° × 0.25° (https://www.earthdata.nasa.gov, last access: 6 September 2024) for the time period of 2017–2023. The retrieved columnar NO₂ is cloud-screened, and the average of the global NO₂ during 2017–2023 was obtained to get an overview of the regions with high NO₂ based on OMI satellite data global measurements as presented in Sect. 2.2.1. These datasets are referred to as OMId (OMI daily) throughout the paper.

2.2 Methodology

2.2.1 Study locations

Taking into account the PGN stations around the globe, and having data availability as specified in Sect. 2.1.2 (version and retrieval code), we selected the co-located AERONET stations with matching latitude and longitude. For multiple co-located AERONET stations, the station having the closest match with the PGN station latitude and longitude, a continuous data flow, and/or larger data availability was selected. By applying these criteria, we identified a total of 33 co-located globally distributed stations to be used for the analysis (Table 1; refer to Table A1 for details regarding station names used by AERONET and PGN and instrument number). These include 11 stations in Europe, 14 in North America and South America, 7 in Asia, and 1 in the Middle East (Fig. 1). Out of these, one station is in the Southern Hemisphere (COM), one is a polar station (NYA), and five are high-altitude (> 1000 m above sea level) stations. Figure 1 also reports the OMId-satellite-based (as described in Sect. 2.1.3) long-term mean of daily NO₂ values between 2017 and 2023, and this shows that the selected stations cover the NO₂ daily mean load representative of conditions ranging from clean (e.g. < 0.2 × 10⁻⁴ mol m⁻²) to polluted (e.g. > 1 × 10⁻⁴ mol m⁻²).

Table 1 Description of the 33 co-located AERONET and PGN stations. The distance of PGN site from AERONET site is mentioned in brackets with sign.

^∗ These sites are co-located (i.e. instruments are in the same building), but the coordinates (latitude, longitude, and altitude) provided in AERONET/PGN have some errors. This is verified with the station Principal Investigators.

Figure 1 (a) Overview of the co-located AERONET and PGN stations and 7-year-averaged (2017–2023) NO₂ (mol m⁻²) from OMId satellite measurements. Panels (b–d) are the focused maps for the clustered locations in North America, Europe, and northeastern Asia, respectively. Sites labelled in red (24 sites) and blue (9 sites) are categorized as urban and rural sites, respectively.

The co-located AERONET and PGN stations have the latitudes of all PGN stations within AERONET latitude ±0.09° and, in most cases, have exactly the same latitudes (Table 1), while the longitudes of the PGN stations are within AERONET's longitude ±0.07° (Table 1). Corresponding to every measurement of AERONET (time of measurement) within a day, the nearest matching PGN measurement (similar time of measurement) was selected and then the PGN data were time-interpolated to the AERONET time stamp for that day. Following this process, we obtained specific comparison data points for each station during the comparison period of 2017–2023, based on the coincident data availability from AERONET and PGN, which are provided in Table 1 (last column). We have categorized all of these stations as either urban or rural sites, based on a simplified assumption that “rural” corresponds to small cities that are in the countryside or adjacent to the ocean, while other sites are “urban”.

2.2.2 NO₂ correction for AOD and AE retrievals

The differences in the OMIc NO₂ used by AERONET for the calculation of AOD from PGN NO₂ VCD (mol m⁻²) are calculated as

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{NO}}_{\mathrm{2}}={\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}-{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}},\end{array}$$

where AERONET OMIc NO₂ is converted from DU to SI (International System of Units) for VCD (which is mol m⁻²) (1 DU = 4.4614 × 10⁻⁴ mol m⁻²) for comparability. AOD is calculated from direct-sun measurements by sun photometers (Cimel sun photometers in the case of AERONET), using the Lambert–Beer law (Eq. 2) that presents the atmospheric attenuation of radiation as

$$\begin{array}{}\text{(2)}& I\left(\mathit{\lambda }\right)=I_{\mathrm{0}}\left(\mathit{\lambda }\right)\times e^{-\left(\begin{array}{c}{m}_{\text{Ray}}{\mathit{\tau }}_{\text{Ray}}+m_{\text{aer}}{\mathit{\tau }}_{\text{aer}}+m_{{\mathrm{O}}_{\mathrm{3}}}{\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}\\ +m_{{\mathrm{NO}}_{\mathrm{2}}}{\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}+m_{{\mathrm{CO}}_{\mathrm{2}}}{\mathit{\tau }}_{{\mathrm{CO}}_{\mathrm{2}}}\\ +m_{{\mathrm{CH}}_{\mathrm{4}}}{\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}}+m_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}{\mathit{\tau }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\end{array}\right)},\end{array}$$

where I(λ) and I₀(λ) represent the radiation intensity at the surface and top of the atmosphere, respectively, at a specific wavelength λ. τ is the total optical depth, and m is the total optical air mass. Total optical depth is the aggregation of the optical-depth contributions from Rayleigh scattering by molecules (τ_Ray), gaseous absorption by ozone (${\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}$), ${\mathrm{NO}}_{\mathrm{2}}\left({\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}\right)$, carbon dioxide (${\mathit{\tau }}_{{\mathrm{CO}}_{\mathrm{2}}}$), methane (${\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}}$), and precipitable water vapour (${\mathit{\tau }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$). m_R, $m_{{\mathrm{O}}_{\mathrm{3}}}$, $m_{{\mathrm{NO}}_{\mathrm{2}}}$, $m_{{\mathrm{CO}}_{\mathrm{2}}}$, $m_{{\mathrm{CH}}_{\mathrm{4}}}$, and $m_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ represent their respective optical air masses, and m_aer is the aerosol optical air mass. The optical air masses are a function of sun elevation. Aerosol optical depth (τ_aer) is calculated from total optical depth (τ) by subtracting the optical-depth contributions from Rayleigh scattering by molecules, gaseous absorption, and/or precipitable water vapour, depending on the wavelength. Here, we only discuss the contribution of NO₂ absorption to AOD and the NO₂ optical-depth estimations (Eq. 3) (Cuevas et al., 2019), calculated as

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}_{{\mathrm{NO}}_{\mathrm{2}}}\left(\mathit{\lambda }\right)}{\mathrm{1000}}}\times {\displaystyle \frac{m_{{\mathrm{NO}}_{\mathrm{2}}}}{m_{\mathrm{a}}}}\times {\mathrm{NO}}_{\mathrm{2}},\end{array}$$

where ${\mathit{\sigma }}_{{\mathrm{NO}}_{\mathrm{2}}}$ is the NO₂ absorption coefficient at wavelength (λ) obtained from Gueymard (1995), and the expression for $m_{{\mathrm{NO}}_{\mathrm{2}}}$ is obtained from Gueymard (1995), while m_a is the optical air mass, and NO₂ VCD is in Dobson units. The NO₂ absorption contribution to the NO₂ optical depth is directly proportional to the NO₂ VCD at a specific wavelength and sun elevation. The bias ΔAOD (or Δτ_aer(λ), as shown in Eq. 5) affecting the AERONET AOD (τ_aer,AERONET) calculation at a specific wavelength is produced by the simplified assumption of OMIc NO₂, and the associated optical depth (which is linear to the NO₂ concentration for an instrument at a specific wavelength and solar elevation; see Eq. 3) is evaluated after exploiting the “real” value of columnar NO₂ from the co-located PGN instrumentation, as shown in Eq. (4) (considering that τ_aer is obtained by subtracting ${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}$ from the total optical depth; hence, ${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}$ is added to τ_aer and newly calculated ${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}}}$ is subtracted to obtain the PGN-corrected τ_aer in Eqs. 4 and 5) as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\tau }}_{\text{aer},\text{PGN}}\left(\mathit{\lambda }\right)=& {\mathit{\tau }}_{\text{aer},\text{AERONET}}\left(\mathit{\lambda }\right)+{\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}\left(\mathit{\lambda }\right)\\ & -\left({\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}\left(\mathit{\lambda }\right)\times {\displaystyle \frac{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}}}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}}\right)\\ =& {\mathit{\tau }}_{\text{aer},\text{AERONET}}-{\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}\left(\mathit{\lambda }\right)\\ & \times \left({\displaystyle \frac{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}}}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}}-\mathrm{1}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\tau }}_{\text{aer}}\left(\mathit{\lambda }\right)=& {\mathit{\tau }}_{\text{aer},\text{AERONET}}\left(\mathit{\lambda }\right)-{\mathit{\tau }}_{\text{aer},\text{PGN}}\left(\mathit{\lambda }\right)\\ =& {\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}\left(\mathit{\lambda }\right)\left({\displaystyle \frac{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}}}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}}-\mathrm{1}\right)\\ =& -{\displaystyle \frac{{\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}\left(\mathit{\lambda }\right)}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}}\left(\mathrm{\Delta }{\mathrm{NO}}_{\mathrm{2}}\right),\end{array}\end{array}$$

where τ_aer,PGN, τ_aer,AERONET, and ${\mathit{\tau }}_{{\mathrm{NO}}_{\mathrm{2}},\text{AERONET}}$ represent the PGN NO₂-corrected AOD, original AERONET OMIc NO₂-based AOD, and OMIc NO₂-based AERONET NO₂ optical depth, respectively (the terms used here are summarized in Table 2). Equation (5) represents the difference in the τ_aer(λ) between AERONET τ_aer and PGN-corrected τ_aer, where the expression for τ_aer,PGN(λ) was obtained from Eq. (4) that led to the second equivalence of Eq. (5), and the third equivalence was obtained using Eq. (1). Therefore, the sign of the AOD bias depends on the sign of ΔNO₂, i.e. the ratio between the OMIc and PGN NO₂. It is also worth noting here that the post-deployment calibrations in level 2.0 data will not have an impact on this analysis of the NO₂-induced differences on AOD differences, as we have considered the relation between the NO₂ difference and AOD difference (Eq. 5) (also from Eq. 3; the NO₂ optical depth is related to columnar NO₂ value, and the other terms will be constant for one instrument at a time stamp or solar elevation and wavelength and are not dependent on the calibration). Therefore, we chose to use level 1.5 data, as described in Sect. 2.1.1, in order to have more comparison points for this analysis. Hence, we define two cases here.

• Case 1. OMIc NO₂ underestimation that is ΔNO₂ < 0, or $\frac{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}}}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}$ > 1, leading to a positive AOD bias (Δτ_aer(λ) > 0) or overestimation of AOD by AERONET (OMIc-based AOD) as compared to PGN-corrected AOD.

• Case 2. OMIc NO₂ overestimation that is ΔNO₂ > 0, or $\frac{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{PGN}}}}{{\mathrm{NO}}_{{\mathrm{2}}_{\mathrm{OMIc}}}}$ < 1, leading to a negative AOD bias (Δτ_aer(λ) < 0) or underestimation of AOD by AERONET (OMIc-based AOD) as compared to PGN-corrected AOD.

Table 2 Summary and description of the terms used in the methodology.

^∗ AERONET is for the Aerosol Robotic Network, PGN is for the Pandonia Global Network, and OMI is for the Ozone Monitoring Instrument.

The spectral variability in AOD is represented by the Ångström exponent (AE), which is obtained from the Ångström power law as

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{aer}}\left(\mathit{\lambda }\right)=\mathit{\beta }\cdot {\mathit{\lambda }}^{-\mathit{\alpha }},\text{(7)}& {\displaystyle }& {\displaystyle }\ln {\mathit{\tau }}_{\text{aer}}\left(\mathit{\lambda }\right)=\ln \mathit{\beta }-\mathit{\alpha }\cdot \ln \mathit{\lambda },\end{array}$$

where α and β represent the AE and the turbidity coefficient, respectively. The negative slope of the least squares regression fit from Eq. (7) is used by AERONET to retrieve AE (Eck et al., 1999) with AOD at all the wavelengths within the considered spectral ranges (here we use all three and four wavelengths within the 340–440 and 440–870 nm wavelength pairs, respectively, for the AE estimations) as

$$\begin{array}{}\text{(8)}& {\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j}=-{\displaystyle \frac{N\sum \ln {\mathit{\tau }}_{\text{aer},i}\cdot \ln {\mathit{\lambda }}_i-\sum {\mathit{\tau }}_{\text{aer},i}\cdot \sum {\mathit{\lambda }}_i}{N\sum (\ln {\mathit{\lambda }}_i{)}^{\mathrm{2}}-(\sum \ln {\mathit{\lambda }}_i{)}^{\mathrm{2}}}}.\end{array}$$

${\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{AERONET}}$ is obtained from AERONET-retrieved AE for two wavelength ranges, namely 340–440 and 440–870 nm. ${\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{PGN}}$ is calculated from the PGN-corrected AOD, i.e. τ_aer,PGN(λ) at wavelengths 340, 380, and 440 nm for the spectral range 340–440 nm; from τ_aer,PGN(λ) at wavelengths 440 and 500 nm; and from τ_aer,AERONET(λ) at 675 and 870 nm for the spectral range 440–870 nm. The difference in the AE is obtained as

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j}={\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{AERONET}}-{\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{PGN}},\end{array}$$

where ${\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j}$ represents the AE in the wavelength range λ_i to λ_j (in our case, these wavelength ranges are 340–440 and 440–870 nm). ${\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{AERONET}}$ and ${\mathit{\alpha }}_{{\mathit{\lambda }}_i-{\mathit{\lambda }}_j,\text{PGN}}$ are the AE based on the AERONET AOD and PGN-corrected AOD, respectively.

2.2.3 AOD and AE trend estimation

We also evaluate the linear trends in AERONET AOD and AE retrievals for a time span of about a decade between 2013–2023 to compare them with the mean AOD and AE differences calculated, as described in Eqs. (5) and (9). The available PGN dataset is for a duration that is quite short, considering the calculation of statistically meaningful trends; hence, we have not considered the trend analysis using PGN-corrected AOD and AE.

The linear AOD and AE trends are evaluated using the weighted least-squares-fitting technique (Weatherhead et al., 1998; Zhang and Reid, 2010; Yoon et al., 2012; Logothetis et al., 2021) as

$$\begin{array}{}\text{(10)}& Y_{\mathrm{m}}=\mathit{\mu }+\mathit{\omega }{X}_{\mathrm{m}}+N_{\mathrm{m}}+S_{\mathrm{m}},\end{array}$$

where m represents the index of the month (m = $\mathrm{1},\mathrm{\dots },M$), M is the total number of months, $M/\mathrm{12}$ is the total number of years, Y_m represents the monthly averaged AOD or AE, X_m represents the decimal number of years since the first month of the time series ($m/\mathrm{12}$), μ represents a constant linear fit offset at the beginning of the time series, ω represents the magnitude of the respective trend per year, and N_m is the residual. The seasonality is taken into account by subtracting S_m, which is the seasonal term calculated as the long-term monthly mean value, from Y_m. For the purpose of deriving statistically significant daily mean values of the aerosol properties (AOD and AE), a minimum of 10 observations on a daily basis was ascertained. Additionally, in order to have a qualified monthly mean, the availability of at least 5 d of measurements on a monthly basis was ensured. The datasets that did not meet these criteria were not considered in the calculation of AOD and AE trends.

The statistical significance of the estimated linear trend (ω) is considered following the methodology presented by Weatherhead et al. (1998), which has been commonly applied for trend detection in AOD by numerous previous studies (e.g. Ningombam et al., 2019; Zhang et al., 2018; Alfaro-Contreras et al., 2017; Adesina et al., 2016; Pozzer et al., 2015; Kumar et al., 2015, 2018; Li et al., 2014; Babu et al., 2013; Hsu et al., 2012), by considering N_m, which follows a first-order autoregressive process as

$$\begin{array}{}\text{(11)}& N_{\mathrm{m}}=\mathit{\phi }{N}_{m-\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{m}},\end{array}$$

where φ is the autocorrelation coefficient (lag−1), and ε_m represents the white noise. The standard deviation of the trend is calculated as

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\mathit{\omega }}\approx {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{N}}}{n^{\mathrm{3}/\mathrm{2}}}}\sqrt{{\displaystyle \frac{\mathrm{1}+\mathit{\phi }}{\mathrm{1}-\mathit{\phi }}}},\end{array}$$

where σ_N represents the standard deviation of N_m, and n is the number of years based on the data availability, taking into account the entire period under consideration (i.e. in our case it is a constant value of 11 years). The trends are considered to be significant when the absolute value of $\mathit{\omega }/{\mathit{\sigma }}_{\mathit{\omega }}$ is above 2.

3 Results and discussion

3.1 Differences between AERONET OMI NO₂ climatology and PGN NO₂ measurements and the impact on AOD measurements

As presented in Sect. 2.2.2, we refer to the OMIc NO₂ underestimation (i.e. ΔNO₂ < 0; $\mathrm{PGN}/\mathrm{OMIc}$ NO₂ ratio > 1) and hence AOD overestimation (ΔAOD > 0) as case 1 and the OMIc NO₂ overestimation (i.e. ΔNO₂ > 0; PGN/OMIc NO₂ ratio < 1) leading to AOD underestimation (ΔAOD < 0) as case 2, which we further discuss here.

Overall, we found 16 (∼ 48 % of all stations) stations in the category of case 1, with mean OMIc NO₂ underestimated as compared to PGN and hence AOD overestimation (Fig. 2a), where 13 locations (∼ 81 % of case 1 stations) are urban sites, and 3 locations (∼ 19 % of case 1 stations) are rural sites. Out of these, six urban stations (DHK, MXC, ATH, LPT, HOU, and ROM; ∼ 37 %) had mean NO₂ underestimation greater than 0.5 × 10⁻⁴ mol m⁻². There were at least 1500 instances with mean ΔNO₂ < −1 × 10⁻⁴ mol m⁻² (Table A2 and A3 in the Appendix) that also showed an AOD overestimation equivalent to or above 0.002. For these cases, the corresponding time series of NO₂ values, differences, and the normalized frequency distribution of the differences are presented in Fig. 3a–f. The mean PGN and OMIc values in DHK are 5.59 × 10⁻⁴ and 1.26 × 10⁻⁴ mol m⁻², respectively, which has higher real (PGN) NO₂ levels, reaching close to 30 × 10⁻⁴ mol m⁻², while OMIc NO₂ remains mostly constant and well within 5 × 10⁻⁴ mol m⁻² (Fig. 3a). In ATH, these values are 2.50 × 10⁻⁴ and 1.20 × 10⁻⁴ mol m⁻², respectively, and in MXC, 3.84 × 10⁻⁴ and 2.01 × 10⁻⁴ mol m⁻², respectively. These stations also have relatively higher real NO₂ values reaching close to 20 × 10⁻⁴ mol m⁻², with OMIc NO₂ being mostly constant at ATH and variable at MXC but well within 5 × 10⁻⁴ mol m⁻² for both stations (Fig. 3b and c). The corresponding AOD differences at 380 nm are 0.015 (∼ 1.0 %), 0.005 (∼ 1.8 %), and 0.007 (∼ 1.7 %) (Table A2 and Fig. A1) for DHK, ATH, and MXC, respectively. At 440 nm, these AOD differences are 0.013 (∼ 1 %), 0.004 (∼ 1.8 %), and 0.005 (∼ 1.7 %) for DHK, ATH, and MXC, respectively (Fig. 2a; Table A2 and Fig. 4). The stations LPT and HOU (Fig. 1), with the NO₂ difference of 0.71 × 10⁻⁴ and 0.58 × 10⁻⁴ mol m⁻², respectively, between OMIc and PGN, showed a mean difference in AOD of 0.003 and 0.002 (∼ 1.1 %) at 380 nm and of 0.002 (∼ 1.1 %) at 440 nm. For ROM, ΔNO₂ was found to be −0.60 × 10⁻⁴ mol m⁻², leading to a mean AOD overestimation of 0.002 at 380 and 440 nm by AERONET OMIc as compared to PGN. LPT, HOU, and ROM has comparatively smaller NO₂ values in the time series (reaching close to 10 × 10⁻⁴ mol m⁻², as per Fig. 3d–f) as compared to stations like DHK and MXC which are located in high-NO₂ zones (as per Fig. 1). The effects of NO₂ differences on AOD at 340 and 500 nm are smaller compared to 380 and 440 nm for all the stations.

Figure 2 NO₂ VCD (mol m⁻²) and AOD differences at 340, 380, 440, and 500 nm for all stations with NO₂ (a) underestimation and (b) overestimation. The NO₂ differences are calculated as OMIc − PGN, and the corresponding AOD differences are calculated as original AERONET AOD minus PGN-corrected AOD (as described in Sect. 2.2.2). The average AOD at each wavelength is plotted as $\mathrm{AOD}/\mathrm{100}$.

Figure 3 (a–j) Time series of NO₂ (mol m⁻²) from OMIc and PGN (black and green dots, respectively), and NO₂ differences (OMIc − PGN; orange dots) and normalized frequency distribution of the NO₂ differences. The 10 panels refer to stations with mean NO₂ differences above 0.5 × 10⁻⁴ mol m⁻² and mean AOD differences above 0.002. The numbers in parentheses represent the mean values.

Figure 4 (a–j) AOD as a function of NO₂ VCD (mol m⁻²). (k–t) AOD differences as a function of AOD at 340, 380, 440, and 500 nm for stations with a mean NO₂ offset of more than 0.5 × 10⁻⁴ mol m⁻² and a mean AOD difference offset above 0.002. For NO₂ underestimation cases (k–p), ΔAOD below 0 for 340 and 500 nm and ΔAOD above 0 for 380 and 440 nm represent positive AOD differences. For NO₂ overestimation cases (q–t), ΔAOD below 0 for 340 and 500 nm and ΔAOD above 0 for 380 and 440 nm represent negative AOD differences.

The underestimation of NO₂ by AERONET OMIc compared to PGN values at stations like DHK and MXC is possibly due to higher pollution levels, which the averaged OMIc climatological interpretation of NO₂ fails to depict, leading to differences compared to the climatological means (Giles et al., 2019). A study by Pavel et al. (2021) on a yearly trend analysis of NO₂ for Dhaka showed a statistically significant positive annual slope of 0.47 ± 0.03 ppb yr⁻¹ for the studied period between 2003 and 2019, which represent an increase in NO₂ levels of ∼ 68 % in 2019 from the base year in 2003, and a similar positive trend was observed by Georgoulias et al. (2019) as 0.29 ± 0.02 × 10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ or 0.05 ± 0.00 × 10⁻⁴ $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ between 1996–2017. The same study by Georgoulias et al. (2019) also revealed a statistically significant positive trend of 0.17 ± 0.09 × 10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ or 0.03 ± 0.01 × 10⁻⁴ $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in NO₂ values for Mexico City.

On the other hand, case 2 had 17 (∼ 52 % of all the stations) stations with the mean NO₂ overestimated by the OMIc when compared to PGN, leading to AOD underestimation (Fig. 2b), with 11 stations (∼ 65 % of the case 2 stations) in the urban area and 6 (∼ 35 % of case 2 stations) in the rural area. Out of these stations, the highest OMIc NO₂ overestimation was observed for four (∼ 23 % of the stations in case 2) urban stations, namely BEI, BRW, TSU, and JYC, with mean differences above 0.5 × 10⁻⁴ mol m⁻², and at least 1500 instances with the overestimation above 1 × 10⁻⁴ mol m⁻² (Table A2 and A3). These four stations also showed the AOD underestimation equal to or above 0.002. The associated NO₂ time series of values, differences, and the normalized frequency distribution of the differences can be found in Fig. 3g–j. The average NO₂ values for BEI were found to be 3.06 × 10⁻⁴ and 4.17 × 10⁻⁴ mol m⁻² from PGN (NO₂ values reaching close to 20 × 10⁻⁴ mol m⁻²; Fig. 3g) and OMIc, respectively; 1.31 × 10⁻⁴ and 1.94 × 10⁻⁴ mol m⁻², respectively, for TSU; 1.54 × 10⁻⁴ and 2.16 × 10⁻⁴ mol m⁻², respectively, for BRW; and 1.75 × 10⁻⁴ and 2.36 × 10⁻⁴ mol m⁻², respectively, for JYC. These differences led to a mean overestimation of NO₂ from OMIc as 1.30 × 10⁻⁴ mol m⁻² for BEI and ∼ 0.62 × 10⁻⁴ mol m⁻² for BRW, TSU, and JYC, which led to an AOD underestimation of ∼ 0.005 for BEI and ∼ 0.002 for BRW, TSU, and JYC.

Stations like BEI showed an overestimation of NO₂ by AERONET OMIc as compared to PGN, possibly due to the reduction in pollution levels as a result of the implementation of environmental protection policies in eastern China (van der A et al., 2017) that may have led to a significant trend reversal of tropospheric NO₂ during the last decade, which OMIc is unable to depict as it considers the average values for time period of 2004–2013. Georgoulias et al. (2019) found a decreasing trend of −1.28 ± 0.78 × 10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ or 0.21 ± 0.13 × 10⁻⁴ $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in tropospheric NO₂ from 2011–2018 (with 2011 being the year of trend reversal from a positive to a negative trend). Another study by Xu et al. (2023) on NO₂ trend analysis in Beijing–Tianjin–Hebei between 2014 and 2020 also revealed a decreasing trend in NO₂ with overall reduction by 44.4 % with reference to the year 2014.

Figure 4 presents the scatterplot of AOD as a function of NO₂ VCD, as well as AOD differences, arising due to NO₂ differences at all considered wavelengths (340, 380, 440, and 500 nm). It is observed that AOD is not correlated with the NO₂ VCD magnitude as observed in Fig. 4a–j, and the AOD differences are also not correlated with the AOD values (Fig. 4k–t). The NO₂ differences are related to the AOD differences, and vice versa, and are not related to the magnitude of AOD or the magnitude of NO₂ VCD as evidenced from Eq. (5).

3.2 Assessment of AOD differences in extreme NO₂ load cases

In this section, we present (Table 3) the scenarios with extreme NO₂ situations; i.e. 10 % highest of the difference cases (from all of the differences presented in Sect. 3.1) are taken into account as percentiles of NO₂ differences with 10 % and 90 % confidence levels for case 1 (NO₂ underestimation by OMIc) and case 2 (NO₂ overestimation by OMIc), respectively (hereafter referred to as the “Extreme” case). Figure 5 presents a comparison of the NO₂ and AOD differences between the Extreme case and whole dataset (hereafter referred to as “All”). It is observed (from Figs. 2 and 5) that the most affected wavelength due to differences in NO₂ absorption representation in AOD calculations is 380 nm, followed by 440, 340, and 500 nm, respectively.

Table 3 Statistics for Extreme cases with 10 % highest NO₂ differences (mol m⁻²) (percentiles (P) at 10 % and 90 % confidence level for case 1 and case 2, respectively).

Figure 5 Comparison of NO₂ VCD (mol m⁻²) and AOD differences (OMIc − PGN) at 340, 380, 440, and 500 nm in Extreme cases, with 10 % of highest NO₂ (a) underestimations and (b) overestimations by OMIc compared to All datasets. The average AOD in the Extreme case at each wavelength is plotted as $\mathrm{AOD}/\mathrm{100}$.

Figure 5a presents the results for case 1, in which the mean differences in an Extreme case were found to be higher than All data case for NO₂ by at least 1 × 10⁻⁴ mol m⁻² and 0.003 for AOD for all stations except NYA and DLG. For the six selected stations from case 1, as discussed in Sect. 3.1, this difference between the Extreme and All case scenario for NO₂ varied from ∼ 2 × 10⁻⁴ mol m⁻², reaching up to even 6 × 10⁻⁴ mol m⁻² (for DHK). The increase in AOD differences for these six stations was found to be above 0.007, reaching up to 0.023 and 0.015 for DHK and MXC, respectively. Similarly, ALD showed ∼ 7 times and ∼ 8 times increase in the differences in NO₂ and AOD, respectively, in the Extreme scenario compared to All datasets.

For case 2, as presented in Fig. 5b, nine stations showed the mean difference between OMIc and PGN NO₂ above 1 × 10⁻⁴ mol m⁻², and the differences in OMIc and the PGN NO₂ difference in the Extreme case from the respective differences in the All dataset were found to reach up to ∼ 2 × 10⁻⁴ mol m⁻². These NO₂ differences led to an average AOD underestimation equivalent to or above 0.002 at 380 and 440 nm at 14 (out of 17) stations by AERONET. The noticeable stations in this case are BEI, JYC, and MNH (Fig. 5b), with the difference in the OMIc and PGN NO₂ difference in the Extreme case from the respective differences in the All dataset being above 1 × 10⁻⁴ mol m⁻², leading to higher AOD differences in the Extreme case than the All dataset by a factor of 0.004 and 0.003 at 380 and 440 nm, respectively. It is noteworthy that for BEI, the mean AOD underestimation between OMIc and PGN reached 0.013 and 0.011 at 380 and 440 nm, respectively, for the mean AOD values of 0.083 and 0.076, respectively. This indicates that high-NO₂ differences in BEI are observed for low AOD cases (Tables 3 and A4), where OMIc overpredicts NO₂ values measured by PGN (Fig. 3g) (Beijing is case 2 of this analysis). Hence, the highest NO₂ differences occur for a low-pollution scenario (i.e. PGN measured NO₂ is lower than OMIc NO₂) and, hence, probably lead to low mean AOD. These cases are about 10 % of what we have considered for Extreme scenario cases for which we have considered the top 10 % of the highest NO₂ differences (for case 1, 90th percentile, and case 2, 10th percentile). Another station to note here is SOL, which showed an increase in the average difference in NO₂, AOD380, and AOD440 from 0.05 × 10⁻⁴ mol m⁻², 0.000 and 0.000 in All datasets (Fig. 5b) to 2.28 × 10⁻⁴ mol m⁻², −0.008 and −0.007, respectively, in the Extreme scenario.

Figure 6 presents the stations with high variations (AOD differences in AERONET from a PGN equivalent to or above 0.005), the mean NO₂ and AOD differences at these stations, and the normalized frequency distribution of the AOD at 340, 380, 440, and 500 nm. A clear shift in the frequency distribution (Fig. 6f1–k4) is observed for Extreme cases – moving away from the All dataset case at the four wavelengths – with a larger shift noticeable at DHK and MXC and a shift in the opposite direction in the case of BEI, which is consistent with the analysis presented in Fig. 5 and Table 3.

Figure 6 ΔNO₂ (mol m⁻²) (a) and ΔAOD at 340, 380, 440, and 500 nm (b–e). The  normalized frequency distribution of AOD differences in the Extreme NO₂ scenario from the whole dataset (referred to as All) for the stations with high variations at corresponding wavelengths (f1–k4).

Figure 7 presents a sensitivity analysis of AOD differences between AERONET and PGN at 380 and 440 nm for all stations, with PGN NO₂ varying between 2 × 10⁻⁴ and 8 × 10⁻⁴ mol m⁻². The median AOD differences are found to be within ±0.01 and go above 0.01 – and even above 0.02 – with the increase in the NO₂ threshold (lower limit) from 2 × 10⁻⁴ to 8 × 10⁻⁴ mol m⁻². Hence, in the case of high-NO₂ loadings, the AOD is expected to have higher uncertainties due to inaccurate NO₂ optical depth estimations.

Figure 7 Variation in AOD differences (AERONET OMIc-based AOD minus PGN-corrected AOD) at 380 and 440 nm for PGN NO₂ varying from (a–d) 2 × 10⁻⁴ to 8 × 10⁻⁴ mol m⁻², respectively, for all stations.

3.3 Effect of climatological vs. real NO₂ values on the Ångström exponent

Due to a differential impact of the NO₂ correction on the spectral AOD, discrepancies between an assumed climatological NO₂ value (OMIc by AERONET) and the real one (PGN based) also impact the AERONET AOD-based computation of the AE. In this section, we present a discussion regarding the differences in the AERONET AOD-based AE and the AE computed from the PGN-corrected AOD, as described in Sect. 2.2.2.

Figure 8 presents the normalized frequency distribution of these AE differences at the wavelength ranges of 340–440 and 440–870 nm. The median of the AE 440–870 nm difference is found to be −0.07 and −0.05 for BEI and BRW, respectively, and within ±0.03 for other stations. The median of the AE 340–440 nm difference is 0.07 for BEI, 0.04 for BRW, and within ±0.03 for the remaining stations. The narrower frequency distribution for stations like DHK can be attributed to the broader AOD distribution (Wagner and Silva, 2008), as shown in Fig. 6f1–f4, and a broader AE distribution at stations like ATH, LPT, HOU, and ROM can be attributed to the narrower AOD distributions at these locations (some examples of AOD distributions are presented in Fig. 6).

Figure 8 Normalized frequency distributions of (a–j) the difference in AE at 440–870 and 340–440 nm retrieved from the AODs based on AERONET OMIc and PGN NO₂. The shaded background area represents the NO₂ underestimation (grey) (a–f) and overestimation (yellow) (g–j) cases.

In AE retrieval, if the AOD relative errors are equal at both wavelengths, then the AE distribution peak reflects the true value, or else there will be a shift in the peak of the AE distribution (Wagner and Silva, 2008). In our case, there is no error at a higher wavelength (870 and 675 nm, as these wavelengths are not affected by NO₂ absorption, and hence PGN NO₂ corrections are not made), and the higher relative positive error at a shorter wavelength (440 and 500 nm) leads to a shift in the peak of the AE difference (ΔAE440–870) distribution towards a positive value. The peak of the distribution of ΔAE340–440 is towards the opposite direction to that of ΔAE440–870 as the error in this case is higher at a higher wavelength (440 nm) than at a lower wavelength (340 nm) in case 1, and a similar but opposite behaviour is observed for case 2. It is also noteworthy that the uncertainty in AE is not very simple to interpret as it is a derivative quantity, and its sensitivity is dependent on the AOD value, as well as any spectral correlations in the AOD uncertainty (Wagner and Silva, 2008; Sayer, 2020). Figure 9 shows the variation in AE differences with NO₂ VCD and AOD values. For NO₂ underestimation cases and with reference to NO₂ VCD (Fig. 9a–f), there is a strong positive bias in AE 440–870 nm (i.e. a higher AE estimation from AERONET compared to a PGN-corrected AOD-based AE estimation) and a negative bias in AE 340–440 nm, while for NO₂ overestimation cases (Fig. 9g–j), the positive and negative biases are not that strongly present as is in the case of NO₂ underestimation. Looking into the AE difference variations with respect to AOD, it was found that high AE differences are associated with low AOD instances.

Figure 9 Scatterplot of the Ångström exponent (AE) difference at 440–870 and 380–500 nm calculated from the AODs based on AERONET OMIc and PGN NO₂-corrected AOD as a function of (a–j) PGN NO₂ VCD (mol m⁻²), and (k–t) AOD at 440 and 380 nm, respectively. The shaded background area represents NO₂ underestimation (grey) and overestimation (yellow) cases.

3.4 Assessment of NO₂ correction on AOD measurements and AE retrievals in rural sites

For the rural sites considered in this analysis, as presented in Figs. 2 and 5, the mean NO₂ underestimation (case 1 as described in Sect. 2.2.2) and overestimation (case 2) between OMIc and PGN was found to be below 0.50 × 10⁻⁴ and 0.40 × 10⁻⁴ mol m⁻², respectively, reaching an underestimation of 1.56 × 10⁻⁴ mol m⁻² for INN and an overestimation of more than 0.40 × 10⁻⁴ mol m⁻² but below 1.00 × 10⁻⁴ mol m⁻² for WAL, BOU, and LDB in the Extreme NO₂ loading scenario. The corresponding impact on AOD mean in case 1 and case 2 was found to be an overestimation and underestimation below 0.002 and 0.001, respectively, at 380 nm and below 0.001 at other wavelengths. Under Extreme NO₂ scenarios, the overestimation reached 0.005 at 380 and 440 nm and 0.004 at 340 nm for INN, while the underestimation was above 0.001 but less than 0.003 for WAL, BOU, and LDB at 380, 440, and 340 nm. The mean AE 440–870 nm difference was found to be positive and within 0.07 for case 1 and negative and within 0.12 for case 2, whereas the mean AE 340–440 nm difference was found to be negative and within 0.06 for case 1 and positive and within 0.07 for case 2.

3.5 Impact of AOD differences on trend analysis

Another aspect of interest relates to the trends in AOD and AE values observed in the last decade, with different magnitudes (and even sign i.e. both overestimation and underestimation cases presented in Sect. 3.1) in different areas of the globe. Hence, in this section, we present the trends based on original AERONET AOD values for a time duration of 2013–2023. In particular, the AOD trends have been calculated based on the AERONET AOD at 380 and 440 nm for stations with larger AOD differences (ΔAOD > 0.002) for the time period between 2013 and 2023 when only considering sites with a data availability of more than 5 years (complete, i.e. all seasons are homogeneously sampled) over this time span.

Table 4 presents the trend analysis using the AERONET AOD and AE. The trends are compared with the mean ΔAOD, which was previously presented in Sect. 3.1. We found two stations with statistically significant negative trends (BEI and JYC) and one with a statistically significant positive trend (HOU) in AOD and negative trends in AE 440–870 nm. HOU, having a positive AOD trend of 0.003 (Table 4), had a mean AOD overestimation of 0.002 at 380 and 440 nm (Table A2), which might have impact on the trends when calculated with the corrected AOD values. Furthermore, the other two stations (BEI and JYC), showing a negative trend in AOD, showed a mean underestimation of AOD as per the analysis presented in Sect. 3.1. It is indicative of how the NO₂ correction could potentially affect realistic AOD trends. The remaining stations (DHK, MXC, ATH, and ROM) could not present a statistically significant trends and, hence, are not discussed here. This analysis signifies the importance of having correct (real) NO₂ values for optical depth calculations that can impact the trend analysis of AOD and AE; however, the true scenario can be unveiled when the trends are calculated with NO₂-corrected AOD.

Table 4 AERONET AOD trend analysis from 2013–2023 at 380 and 440 nm.

3.6 Pandora NO₂ vertical column density spatial representativeness

In this section, we try to look into the spatial representativeness of the Pandora instruments for the locations, as discussed in previous sections. Figure 10 shows the 7-year-averaged OMId satellite values based on the spatial distribution of NO₂ VCD (also presented in Fig. 1), and the statistics are presented in Table 5. The location of Pandora (marked with red dots) represents the centre of the circular area (red circles), which is considered according to the OMI satellite overpass (yellow dots). The differences are calculated based on the area-averaged NO₂ values from OMId satellite and PGN measurement averages. For stations like DHK and MXC that have higher-NO₂ values, the area-averaged differences increase with the increase in the area, while other stations like ATH, LPT, HOU, and ROM showed a comparatively smaller variation in the differences. For BEI, the differences were constant until second circular area around the Pandora site and then increased with the increasing radius and showed a maximum difference for the outermost circle.

Figure 10 Spatial variation in NO₂ VCD from OMI (7-year-averaged value, as presented in Fig. 1, i.e. during 2017–2023). The red (at the centre) and yellow dots represent the PGN location and the satellite overpass, respectively. The red circles centred around the PGN location are calculated with a radius representative of the distance between the PGN location and satellite overpass.

Table 5 Average NO₂ VCD PGN − OMId satellite difference (× 10⁻⁴ mol m⁻²) of circles centred at the PGN site and of the radius increasing as per the difference between PGN site and OMI satellite overpass. The circles represent the area around the centre and are numbered according to the increasing distance from the centre. The values in parentheses represent the difference in the average NO₂ values of the respective circle from circle 1.

For sites with homogeneous NO₂ distributions, a Pandora instrument can be considered for VCD for a larger surrounding area, while for the regions with less homogeneous NO₂ distributions, there can be a limited representation of NO₂ in the surrounding area by a Pandora instrument (Liu et al., 2024). Moreover, closely located PGN sites like LPT and HOU can be used to include the regional spatial variation in the NO₂. In our analysis, these two closely located stations of LPT and HOU (Fig. 1), with the NO₂ difference of 0.71 × 10⁻⁴ and 0.58 × 10⁻⁴ mol m⁻², respectively, between OMIc and PGN, showed a mean difference in AOD of 0.003 and 0.002 (∼ 1.1 %) at 380 nm, respectively, and of 0.002 (∼ 1.1 %) at 440 nm. Also, Drosoglou et al. (2023a) analysed the spatiotemporal variability in NO₂ for ATH by synergistically using Pandora and satellite (TROPOMI – TROPOspheric Monitoring Instrument) observations. Thus, this could be another approach, such as using a high-resolution satellite VCD for NO₂ characterization for real-time NO₂ estimations or for the improvement of the climatology used for NO₂ optical depth estimation.

4 Conclusion

This work was based on the Drosoglou et al. (2023b) findings showing the NO₂ effects on AOD measurements for Rome, Italy. Here we tried to expand the investigation to all stations with co-located PGN Pandora and AERONET Cimel instruments. We present the analysis of NO₂ differences between AERONET OMI climatology and the PGN dataset, focusing on the assessment of the impact on AOD at 340, 380, 440, and 500 nm from 33 worldwide co-located AERONET and PGN stations. About half of these stations (∼ 81 % of which are in the urban area and the remainder in the rural area) showed an underestimation of NO₂ values by AERONET OMI climatology compared to the real (PGN) NO₂ measurements that could be possibly due to higher pollution levels that averaged AERONET OMI climatological interpretations of NO₂ fail to depict. While the other stations (∼ 65 % of which were urban sites and the remainder were rural sites) showed an overestimation of NO₂ which could possibly be due to the reduction in pollution levels as an outcome of the implementation of environmental protection policies (in the last decade) that may have led to a significant NO₂ trend reversal which AERONET OMI climatology might not be able to depict due to the fact that it considers the average values for time period of 2004–2013.

The correction in AERONET AOD based on PGN NO₂ showed differences from the AERONET OMI climatology-based AOD. The analysis was further focused on 10 stations that showed a minimum mean NO₂ and AOD (at 380 and 440 nm) differences of 0.5 × 10⁻⁴ mol m⁻² and 0.002, respectively. Among these, six stations (DHK, MXC, ATH, LPT, HOU, and ROM) belonged to case 1, having the underestimation of NO₂ and overestimation of AOD, while four stations (BEI, TSU, BRW, and JYC) showed the overestimation of NO₂ leading to AOD underestimation (case 2). The AOD bias was found to be the most affected at 380 nm due to NO₂ differences followed by 440, 340, and 500 nm, respectively.

Further assessment of AOD differences in Extreme NO₂ loading scenarios (i.e. 10 % highest difference instances taken into account as percentiles of NO₂ differences with 10 % and 90 % confidence levels for case 1 and case 2, respectively) revealed higher AOD differences in all cases, with a much more significant increase in the 10 stations mentioned above, along with 3 more stations (ALD, SOL, and MNH), compared to their respective All dataset mean AOD differences. Furthermore, the sensitivity analysis based on the PGN NO₂ variation from 2 × 10⁻⁴ to 8 × 10⁻⁴ mol m⁻² revealed that in case of high-NO₂ loadings, the AOD is expected to have higher uncertainties due to an inaccurate NO₂ optical depth representation by AERONET OMI climatology.

Due to the impact of the NO₂ correction (discrepancies between the AERONET OMI climatological representation of NO₂ values and the real NO₂ measurement values by PGN) on the spectral AOD, the AOD-derivative product, AE, is also impacted. The normalized frequency distribution of AE was found to be narrower for a broader AOD distribution for some stations, and vice versa for other stations. For the wavelength pair used in AE estimation, a higher relative AOD error at the shorter wavelength led to a shift in the peak of the AE distribution towards a positive value, and a higher relative AOD error at a higher wavelength led to the shift in the peak of the AE distribution towards a negative value for the AOD overestimation case, and vice versa for the AOD underestimation case. Also, it is to be noted that the uncertainty in the AE is difficult to interpret due to the AE being a derivative quantity, and its sensitivity depends on the AOD value, as well as any spectral correlations in the AOD uncertainty.

The rural locations considered in this analysis showed mean NO₂ differences mostly below 0.50 × 10⁻⁴ mol m⁻² for case 1 and case 2. AOD differences were found to be mostly below 0.001 at all wavelengths except 380 nm, which had these differences below 0.002. Slightly higher (compared to the All dataset scenario for rural locations) NO₂ and AOD differences were observed in Extreme NO₂ loading scenarios to about 1.50 × 10⁻⁴ mol m⁻² and 0.005, respectively, for some stations.

An AOD and AE trend assessment was made for about a decade for stations with AOD differences above 0.002 and with more than 5 years of data availability, based on the original (i.e. AERONET OMI climatological NO₂) AERONET AOD. Stations having comparable mean AOD overestimations or underestimations with the estimated trends revealed that if the trends can be calculated for these stations with the NO₂-corrected AOD, then there can be impacts on the trend values. This analysis is an indication of how the NO₂ correction could potentially affect realistic AOD trends. However, the true scenario can be unveiled only with the trends that are calculated with NO₂-corrected AOD values. For future analysis, it would be interesting to see how the NO₂-based AOD correction would impact the AOD and AE trends; i.e. how much would the trends deviate when using the corrected AODs?

In general, average AOD-related overestimation or underestimation due to differences in the actual and climatological NO₂ inputs are low, with the exception of few stations for which a decade-old satellite-based NO₂ climatology fails to capture the local NO₂ variability and its absolute levels. However, in the case of high-NO₂ events (days), such differences are important, as for the top 10 % of high-NO₂ cases (these high-NO₂ difference cases are not associated with high AOD cases but are related to high levels of pollution and/or changes in the pollution trends in the past decade; Fig. A1 and Table A4); for 10 of the stations, the impact on the AODs is close to the limit or higher than the reported 0.01 uncertainty reported by Giles et al. (2019) and Eck et al. (1999) for AERONET AOD measurements. Taking into account that this uncertainty is a result of various aspects such as calibration (primarily), postprocessing, and instrument/measurement uncertainty, the NO₂-related contribution can be considered relatively significant. Higher spatial and temporal resolution and updated NO₂-satellite-based climatology or the use of co-located Cimel–Pandora retrievals could limit the reported NO₂-related AOD uncertainties, especially in urban areas where NO₂ can be highly variable.

Moreover, some AOD measuring networks (e.g. SKYNET, Nakajima et al., 2020; GAW-PFR, Kazadzis et al., 2018) do not officially take into account the NO₂ optical depth in AOD measurements, and in this case the NO₂ correction will be considered a systematic overestimation of AOD. For the GAW-PFR network, the NO₂-absorption-based error in AOD measurements can be assumed to be negligible as the GAW remote stations have low-NO₂ concentrations (the annual mean values of NO₂ optical depth are < 0.001 in general; Kazadzis et al., 2018). However, it might be of some significance for stations located in polluted areas, especially in Asia, or during extreme events such as wildfires, which are becoming more frequent as a consequence of climate change. As a future endeavour, it would also be interesting to look into the impact of NO₂-based corrections on AOD and other aerosol property retrievals, especially in ground-based aerosol remote sensing stations located in high-pollution zones such as those of SKYNET, which has established regional sub-network groups in China, Europe, India, Japan, South Korea, Mongolia, and Southeast Asia. Finally, technological improvements and the wide range of instrumentations such as real-time NO₂ monitoring from the Pandonia global network, high-spatial-resolution real-time satellite-based observations (such as TROPOMI), and the foreseen high-temporal-resolution NO₂ products (such as from Sentinel 4 and TEMPO satellites) could be directly used for contributing towards the improvement of aerosol property retrievals, specifically in the spectral range (∼ 340–500 nm), which is significantly affected by NO₂ absorption.

This analysis highlights the importance of accurate NO₂ optical depth representation with the best possible scenario (i.e. high-frequency and accurate available NO₂ measurements from Pandora instruments); however, concerning the implementation into the global AOD networks (such as AERONET, GAW-PFR, or SKYNET), the utilization of satellite data is required to account for all the stations in the network.

Appendix A

List of acronyms​​​​​​​

τ	Optical depth
α	Ångström exponent
λ	Wavelength
Δ	Difference
AE	Ångström exponent
AERONET	Aerosol Robotic Network
AOD	Aerosol optical depth
DU	Dobson unit
GAW-PFR	Global Atmospheric Watch – Precision Filter Radiometer Network
OMI	Ozone Monitoring Instrument
OMIc	OMI climatology
OMId	OMI daily
PGN	Pandonia Global Network
TEMPO	Tropospheric Emissions: Monitoring of Pollution
TROPOMI	TROPOspheric Monitoring Instrument
VCD	Vertical column density

Table A1 AERONET and PGN co-located station information.

^∗ These sites are co-located (i.e. instruments are in the same building) but the coordinates (latitude, longitude, and altitude) provided in AERONET/PGN have some errors at the time of submission. This is verified with the stations' principal investigators.

Table A2 NO₂ (mol m⁻²), AOD (380 and 440 nm), and AE (440–870 nm) differences. All differences are as OMIc − PGN.

Table A3 NO₂ (mol m⁻²), AOD (340 and 500 nm), and AE (340–440) differences. All differences are as OMIc − PGN.

Figure A1 AERONET (a) AOD and (b) AOD percentage difference as a function of NO₂ VCD for 10 % highest NO₂ cases for 10 stations (DHK, MXC, ATH, LPT, HOU, ROM, BEI, TSU, BRW, and JYC).

Table A4 Comparison between NO₂ optical-depth-based bias and relative percentage differences in AOD at 380 nm in Extreme NO₂ cases.