Horizontal distribution of tropospheric NO₂ and aerosols derived by dual-scan multi-wavelength multi-axis differential optical absorption spectroscopy (MAX-DOAS) measurements in Uccle, Belgium

Dual-scan ground-based multi-axis differential optical absorption spectroscopy (MAX-DOAS) measurements of tropospheric nitrogen dioxide (NO₂) and aerosols were carried out in Uccle (50.8^∘ N, 4.35^∘ E; Brussels region, Belgium) for 2 years from March 2018 to February 2020. The MAX-DOAS instrument operated in both UV and visible wavelength ranges in a dual-scan configuration consisting of two submodes: (1) an elevation scan in a fixed viewing azimuthal direction and (2) an azimuthal scan in a fixed low elevation angle (2^∘). By analyzing the O₄ and NO₂ differential slant column density (dSCD) at six different wavelength intervals along every azimuthal direction and by applying a new optimal-estimation-based inversion approach (the so-called mapping MAX-DOAS technique), the horizontal distribution of the NO₂ near-surface concentrations and vertical column densities (VCDs) as well as the aerosol near-surface extinction coefficients are retrieved along 10 azimuthal directions. The retrieved horizontal NO₂ concentration profiles allow the identification of the main NO₂ hotspots in the Brussels area. Correlative comparisons of the retrieved horizontal NO₂ distribution were conducted with airborne, mobile, air quality model, and satellite datasets, and overall good agreement is found. The comparison with TROPOMI observations from operational and scientific data products reveals that the characterization of the horizontal distribution of tropospheric NO₂ VCDs by ground-based measurements and an adequate a priori NO₂ profile shape in TROPOMI retrievals lead to better consistency between satellite and ground-based datasets.

1 Introduction

Aerosols and nitrogen dioxide (NO₂) play a crucial role in tropospheric chemistry. NO₂ is an important tropospheric pollutant mainly emitted by combustion processes and nitrogen fertilizers used in agriculture (Seinfeld and Pandis, 2016). Traffic, domestic heating, industrial activities, and power plants are the largest NO₂ emitters (Tack et al., 2021). Beyond its harmful effects on human health (Chen et al., 2007), NO₂ participates in the formation of tropospheric ozone (O₃) by a nonlinear photochemical mechanism, which involves volatile organic compounds (VOCs).

Aerosols with a small diameter can penetrate deeply into the lungs, causing millions of premature deaths around the world per year (Khomenko et al., 2021). Additionally, aerosols influence the Earth's climate system by changing its radiation budget by scattering and absorbing sunlight (Quaas et al., 2008). In the boundary layer of urban regions, the horizontal distribution of NO₂ is highly heterogeneous given the fact that it is a short-lived species (Beirle et al., 2003). For those reasons, the regional and global monitoring of NO₂ and aerosols at high spatial resolution is crucial.

Since 1995, with the ERS-2 GOME (Global Ozone Monitoring Experiment) instrument (Burrows et al., 1999), satellite nadir air quality measurements of atmospheric backscattered sunlight in the UV–visible range have provided daily global tropospheric column measurements of numerous trace gases, such as NO₂. Many satellite missions dedicated to air quality monitoring followed over the next years with increasing spatial resolution. More recently, the TROPOspheric Monitoring Instrument (TROPOMI) launched on board the Sentinel-5P Precursor (S5P) platform in October 2017 reached an initial spatial resolution of 7×3.5 km² and was augmented on 6 August 2019 to 5.5×3.5 km². Due to TROPOMI's fine spatial resolution, monitoring the horizontal distribution of NO₂ in urban regions and identifying specific emission sources are made easier than with previous satellite missions, but TROPOMI still cannot fully capture the fine-scale (sub-kilometer) structures in the effective NO₂ field. Consequently, TROPOMI requires further attention concerning its measurement validation.

Tropospheric vertical columns of many trace gases like NO₂, formaldehyde (HCHO), sulfur dioxide (SO₂), nitrous acid (HONO), and O₃ can be retrieved by the multi-axis differential optical absorption spectroscopy (MAX-DOAS) technique (Hönninger et al., 2004; Wittrock et al., 2004; Pinardi et al., 2008, 2013; Clémer et al., 2010; Hendrick et al., 2014; Irie et al., 2011, 2012; Sinreich et al., 2007; Wagner et al., 2011; Wang et al., 2018). In recent years, MAX-DOAS measurements have been widely used as reference datasets for the validation of nadir airborne and spaceborne air quality measurements. MAX-DOAS instruments measure the scattered sunlight in the UV and visible spectral ranges at multiple elevation angles above the horizon. For absorbers located close to the surface, such as tropospheric NO₂, higher sensitivity is achieved for low MAX-DOAS elevation angles. During the last years, MAX-DOAS measurements in more than one azimuthal direction have emerged (Ortega et al., 2016; Wang et al., 2014; Chan et al., 2020; Schreier et al., 2021). Multi-azimuthal MAX-DOAS measurements offer many possibilities for air quality monitoring, such as a better characterization of the effective NO₂ field around the station. These ground-based datasets can be valuable for validating satellite missions with fine spatial resolution in regions where the NO₂ horizontal distribution is heterogeneous, such as urban and suburban areas.

In this study, a new aerosol and NO₂ horizontal distribution inversion approach based on 2 years (March 2018–February 2020) of dual-scan multi-wavelength MAX-DOAS measurements in Uccle (Brussels-Capital Region, Belgium) is presented. We refer to this new approach as the mapping MAX-DOAS technique. In every azimuthal viewing direction, parameterized NO₂ near-surface concentrations, NO₂ tropospheric columns, and aerosol extinctions measured at six different wavelengths are used as input in a new horizontal distribution inversion approach. On this basis, the near-surface aerosol extinction and NO₂ horizontal distributions are retrieved at a spatial resolution of about 3 km in a range of about 20 km around the measurement site. These horizontal profiles are used to validate co-located TROPOMI tropospheric NO₂ columns. We use (i) the operational TROPOMI data product, (ii) a so-called diagnostic dataset, which is a small sample but with major updates compared to the operational product (van Geffen et al., 2022), and (iii) a European TROPOMI product for which profiles from the CAMS (Copernicus Atmospheric Monitoring Service) regional chemistry transport model (CTM) ensemble (S5P-CAMS) are used (Douros et al., 2022). One complete year of data (March 2018–March 2019) and two wavelength intervals (one in the UV and one in the visible) have already been used in Dimitropoulou et al. (2020). It is proven that multi-azimuthal (the so-called dual-scan) MAX-DOAS measurements significantly improve the agreement between ground-based and TROPOMI tropospheric NO₂ column observations over the Brussels area. By adding the multi-wavelength aspect, the present work represents an extension of the former study.

The paper is organized into seven sections: in Sect. 2, the measurement site with the MAX-DOAS experimental setup and the multi-wavelength DOAS analysis are presented. In Sects. 3 and 4, the TROPOMI tropospheric NO₂ measurements and ancillary measurements used in this study are described, respectively. Section 5 is composed of two main parts: Sect. 5.1 is a detailed description of the dual-scan multi-wavelength MAX-DOAS retrieval method, and Sect. 5.2 is the horizontal aerosol and NO₂ distribution inversion approach (i.e., the mapping MAX-DOAS technique). In Sect. 6, main results followed by correlative comparisons of the retrieved ground-based and satellite horizontal NO₂ distribution are presented. Finally, in Sect. 7, conclusions and future perspectives are given.

2 Dual-scan multi-wavelength MAX-DOAS measurements

2.1 Measurement site and experiment setup

The Brussels-Capital Region is the most densely populated area in Belgium, where pollutant concentrations, such as NO₂, are often high because of anthropogenic activities (Tack et al., 2021).

A MAX-DOAS dual-scan instrument was operated by BIRA-IASB (Koninklijk Belgisch Instituut voor Ruimte-Aeronomie – Institut Royal d'Aeronomie Spatiale de Belgique) in Uccle from January 2017 to February 2020. Uccle is located to the south of the city center of Brussels and to the west of a large forested area (Bois de la Cambre). Therefore, it is an ideal site to perform MAX-DOAS observations under moderate to high pollution conditions. Additionally, the characterization of the horizontal distribution of NO₂ and aerosols at high spatial resolution is of great interest here because of the heterogeneity of the pollution sources (car traffic, national airport, power plant) in the capital region of Brussels.

The MAX-DOAS dual-scan instrument is composed of the following parts: the optical head mounted on a sun tracker, two spectrometers (UV and visible) inside a thermo-regulated box, and the data acquisition unit. The optical head and the two spectrometers are connected with optical fibers. A more detailed description of the BIRA-IASB MAX-DOAS dual-scan instrument can be found in Dimitropoulou et al. (2020).

From March 2018 to February 2020, the MAX-DOAS instrument operated in a dual-scan viewing mode. Two different submodes compose one complete measurement scan (see Dimitropoulou et al., 2020): (1) a vertical scan in nine different elevation angles (EAs) in one fixed telescope azimuth angle (TAA; northeast direction, i.e., towards the city center and the national airport) and (2) a horizontal scan in nine different azimuthal directions at a fixed elevation angle (2^∘ above the horizon).

Several azimuthal viewing directions were tested to obtain an optimal horizontal sampling without any obstacles in the different viewing directions (see Fig. 1). The selection of more azimuthal directions towards the north, northeast, and northwest was made considering the location of the main NO₂ emission sources and, consequently, the highly variable NO₂ horizontal distribution towards these directions.

Figure 1 The experimental setup of the BIRA-IASB dual-scan MAX-DOAS instrument. Each line is color-coded according to the different setups that were used from March 2018 to February 2020. The length of each line is equal to 20 km, which corresponds to the typical horizontal sensitivity for the MAX-DOAS measurements in the present study (see Fig. 18). The black square shows the MAX-DOAS instrument location, the black polygon the national airport, the black dots the NO₂ hotspots emitting more than 10 kg NO_x h⁻¹ (Emission Inventory of the Belgian Interregional Environment Agency, 2017; Frederik Tack, personal communication, 2019), and the black line the Brussels Ring motorway.

The integration time for each measured radiance spectrum is 60 s, resulting in a full scan duration of approximately 20 min.

2.2 Multi-wavelength DOAS analysis

The measured radiance spectra of a full measurement scan are analyzed using the QDOAS spectral fitting software developed by BIRA-IASB (Fayt et al., 2017). The DOAS technique separates the narrow absorption features of trace gases in the UV–visible spectral range from a spectral background caused mainly by Mie and Rayleigh scattering and instrumental effects. The trace gas concentration integrated along the light path in a measured spectrum relative to the amount of the same absorber in a reference spectrum is the primary product of the DOAS analysis and is called differential slant column density (dSCD). Here, average zenith spectra before and after each measurement scan are used as a reference.

The O₄ and NO₂ dSCDs are retrieved in six different wavelength intervals: three intervals in the UV spectral range (330–361, 350–370, and 360–383.5 nm) and three in the visible range (420–460, 450–490, and 510–540.1 nm). These fitting windows were selected to optimize the determination of the O₄ and NO₂ dSCDs at the maximum number of different O₄ absorption bands available in the wavelength domain of the instrument. Figure 2 shows an example of the O₄ and NO₂ fits in all the intervals used in the present work. In each chosen fitting window, we select a reference wavelength, which corresponds to the maximum of an O₄ absorption peak (or close to it) in the respective wavelength intervals (see Fig. 2), and it is subsequently used for radiative transport calculations and further analysis. The different reference wavelengths are 343, 360, 380, 447, 477, and 530 nm (see Fig. 2). To optimize the derivation of the dSCDs at the six selected wavelengths, the fit of a slope parameter, which accounts for the variation of the dSCD within the fitting interval (Puķīte et al., 2010), is necessary. This is especially important when the reference wavelength is not located in the center of the fitting window (i.e., 330–361, 350–370, 420–460, and 450–490 nm). The DOAS settings used for each fitting interval are presented in Table S1 in the Supplement. As shown in this table, two different O₄ cross-sections are used in this study: (1) Finkenzeller and Volkamer (2022) in the UV fitting intervals and (2) Thalman and Volkamer (2013) in the visible fitting intervals. The main motivation for using the O₄ cross-section from Finkenzeller (measured at 25 ^∘C) in the UV fitting intervals is the significant improvement of the fit quality and the reduction of the uncertainties for the UV retrievals. Sensitivity tests and comparisons with radiative transport simulations also show that the resulting O₄ and NO₂ dSCDs are consistent throughout the whole wavelength range covered by the six intervals. For NO₂ and O₃, which are the strongest absorbers in all the fitting windows, a correction for the solar I₀ effect (Aliwell et al., 2002) is applied. A high-resolution solar atlas (Kurucz et al., 1984) is used for the wavelength calibration of the measured spectra.

Figure 2 Results of the O₄ and NO₂ fit for the six selected fitting windows from the dual-scan MAX-DOAS measurements in Uccle (2 June 2019 at 07:05 UTC). The measured spectra are represented with black lines, while the fit results are shown with red lines. The blue lines represent the six reference wavelengths.

3 TROPOMI tropospheric NO₂ measurements

In the present study, MAX-DOAS tropospheric NO₂ VCDs are used to validate co-located TROPOMI satellite observations. TROPOMI is a passive grating imaging spectrometer flying on board the S5P satellite platform. It covers the UV–visible (250–500 nm), near-infrared (710–770 nm), and shortwave infrared (2314–2382 nm) spectral ranges (Veefkind et al., 2011). TROPOMI measures in a push-broom configuration with a full swath width as wide as 2600 km, and it provides daily global coverage at a spatial resolution (true-nadir pixel size) of 7×3.5 km², which is further improved to 5.5×3.5 km² on 6 August 2019. The TROPOMI tropospheric NO₂ algorithm has been developed at KNMI and uses a retrieval assimilation modeling system that is based on the 3-D global TM5 chemistry transport model (van Geffen et al., 2019; Williams et al., 2017).

We use the reprocessed (RPRO) and offline (OFFL) datasets of the TROPOMI L2 tropospheric NO₂ column product (see Table 1 for the corresponding versions). According to the guidelines provided by van Geffen et al. (2019), RPRO datasets are available only for the first period of the present study (see Table 1). For the remaining periods, OFFL datasets are used, which are the main data products available within 2 weeks from the TROPOMI measurement. To ensure the best measurement quality, only pixels with a quality assurance value larger than 0.75 are used. This quality flagging eliminates pixels with a cloud radiance fraction larger than 0.5, snow or ice, and erroneous retrievals (Eskes and Eichmann, 2022).

Table 1 TROPOMI NO₂ processor versions used in the present study.

Next to operational products, two additional TROPOMI datasets are also used (see Sect. S2 in the Supplement). In the first one, the TROPOMI retrieval is performed with different a priori profiles (Douros et al., 2022). The coarse TM5-MP a priori NO₂ profiles, using a spatial resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, are replaced by NO₂ profile shapes from the CAMS (Copernicus Atmospheric Monitoring Service) regional chemistry transport model (CTM) ensemble at a spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ (S5P-CAMS product). The replacement of coarse a priori information with finer information can lead to significant changes in the TROPOMI-retrieved NO₂ tropospheric columns. The available dataset covers October 2018 to March 2020 (OFFL dataset, L2, and version 01.02.00 up to 01.04.00; Eskes and Eichmann, 2022).

In the second additional product, the TROPOMI retrieval is performed with an improved cloud product (Eskes et al., 2022; van Geffen et al., 2019). According to van Geffen et al. (2019), the improvement of the FRESCO-S cloud pressure retrieval scheme to the FRESCO-wide product has an impact on the NO₂ air mass factors (AMFs) and consequently on the NO₂ tropospheric columns over polluted areas. More precisely, the existing FRESCO-S product had a negative bias in the cloud height values, which resulted in a low NO₂ tropospheric column (Compernolle et al., 2021). The TROPOMI tropospheric NO₂ columns are retrieved using an improved FRESCO-S cloud retrieval scheme, called FRESCO-wide, in v1.4 since 29 November 2020 (Eskes and Eichmann, 2022). In the present study, the diagnostic datasets (DDSv2) are used, which are an ensemble of reprocessed data for past periods analyzed with new versions (van Geffen et al., 2022). Over the MAX-DOAS measurement time period, only DDSv2 data corresponding to OFFL datasets (v1.2 and v1.3) are available. Excluding the spin-up period needed by TM5-MP, only four data periods are available for our comparisons (i.e., 30 June–6 July 2018, 30 December 2018–5 January 2019, 30 March–5 April 2019, and 17–23 September 2019).

4 Ancillary measurements

In the present study, ancillary measurements are used for two main purposes: (1) as a priori information on the retrieval of the MAX-DOAS NO₂ and aerosol horizontal profiles as well as (2) to validate the retrieved MAX-DOAS NO₂ horizontal profiles.

First, the CIMEL Sun photometer is a multi-channel photometer, which measured the direct solar irradiance and sky radiance at the Earth's surface in discrete wavelengths in the UV, visible, and near-IR wavelengths of the spectrum. Aerosol optical thickness measurements at the Brussels-Capital Region from a co-located CIMEL Sun photometer of AERONET (https://aeronet.gsfc.nasa.gov/, last access: 11 July 2022) are used in the retrieval of the MAX-DOAS aerosol horizontal profiles (see Sect. 5.3).

Secondly, independent NO₂ horizontal profiles, used as a priori information in the new OEM-based (OEM: optical estimation method) horizontal distribution inversion approach, are provided by the RIO model. RIO is a land use regression model based on the interpolation of the hourly NO₂ near-surface concentrations measured by the in situ telemetric air quality network in Belgium (Hooyberghs et al., 2006; Janssen et al., 2008; see https://irceline.be/en, last access: 13 July 2022). RIO provides hourly NO₂ concentration maps at a 4×4 km² spatial resolution.

To validate the retrieved MAX-DOAS NO₂ horizontal profiles, independent measurements are necessary in the Brussels-Capital Region.

For the S5P validation campaign over Belgium (S5PVAL-BE, https://s5pcampaigns.aeronomie.be/, last access: 11 July 2022), airborne measurements in the two largest urban regions over Belgium, i.e., Antwerp and Brussels, took place from 26 to 29 June 2019 (Tack et al., 2021). The Airborne Prism EXperiment (APEX) imaging spectrometer was used to measure the horizontal distribution of tropospheric NO₂ columns with a spatial resolution of approximately 75 m×120 m (Tack et al., 2017, 2019). During the S5PVAL-BE flight over Brussels, car-mobile DOAS observations were performed by the BIRA-IASB mobile DOAS, the so-called AEROMOBIL (Merlaud, 2013). The AEROMOBIL consists of a compact double Avantes spectrometer simultaneously recording scattered light in two elevation angles (i.e., one at 30^∘ elevation angle and one at zenith). The AEROMOBIL was used to measure the spatial distribution of tropospheric NO₂ columns mainly over the Ring road of Brussels. These airborne and car-mobile measurements are compared with the retrieved MAX-DOAS NO₂ horizontal profiles (see Sect. 6.2).

5 Description of the mapping MAX-DOAS technique

First, the measured radiance spectra in the UV and visible wavelength ranges are analyzed in the six different fitting windows listed in Sect. 2.2. Then, the optimal-estimation-method-based (OEM-based) Mexican MAX-DOAS Fit (MMF) algorithm is applied to the O₄ and NO₂ dSCDs in the main azimuthal direction (and at 477 nm) to retrieve vertical NO₂ profiles and obtain information about the vertical extent of NO₂ in the troposphere via the mixing layer height of NO₂ (MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$; see Sect. 5.1).

As an intermediate step, radiative transfer model (RTM) simulations are performed (see Table 2 and Sect. 5.2) to obtain information about the horizontal sensitivity ($L_{{\mathrm{NO}}_{\mathrm{2}}}$) and aerosol optical depth (AOD) as a function of O₄ dSCDs, wavelength, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$.

Table 2 RTM inputs for the simulations of $L_{{\mathrm{NO}}_{\mathrm{2}}}$ at the six selected wavelengths (343, 360, 380, 447, 477, and 530 nm).

Then, in the next step, a new dual-scan parameterization technique is applied to the O₄ and NO₂ dSCDs at the six different wavelengths and in all the azimuthal directions with MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, measured O₄ dSCDs, and measurement geometry being the main input parameters to retrieve the horizontal sensitivity of NO₂ and, consequently, the NO₂ near-surface concentrations, VCDs, and near-surface aerosol extinction as a function of distance from the instrument (see Sect. 5.2).

In the final step, a new OEM-based horizontal distribution inversion approach is developed using the six near-surface NO₂ concentrations and aerosol extinction values per azimuthal direction to retrieve NO₂ and aerosol extinction horizontal profiles in an output grid of 500 m thickness (see Sect. 5.3).

A flowchart describing the mapping MAX-DOAS technique is shown in Fig. 3.

Figure 3 Mapping MAX-DOAS technique flowchart.

5.1 Aerosol and NO₂ OEM-based profile retrievals

The optimal-estimation-based MMF inversion algorithm (Friedrich et al., 2019) is applied to retrieve the aerosol extinction coefficient and NO₂ vertical profiles for each MAX-DOAS elevation scan in the main azimuthal direction at 360 and 477 nm. First, the O₄ measurements are used to retrieve the aerosol extinction profile. Several studies indicated the importance of applying a scaling factor (≠1) to the observed O₄ dSCDs to bring them in agreement with simulated O₄ dSCDs by radiative transfer modeling (Wagner et al., 2009; Clémer et al., 2010; Merlaud et al., 2011; see also Wagner et al., 2019, Table 1 for a comprehensive list of all those studies). However, there is no consensus on the fundamental reason for applying this scaling (see, e.g., Ortega et al., 2016). As found by Tirpitz et al. (2021), the choice of the scaling factor has only a small effect on the performance of the trace gas retrieval, so we decided not to apply it in the present study. The aerosol extinction profile retrieved from each scan is used as an input to the radiative transfer calculations used to retrieve the NO₂ profile. Further details about the MMF inversion algorithm, the input a priori parameters, the quality check of each scan, and the estimated uncertainties of the aerosol and NO₂ vertical profile can be found in Dimitropoulou et al. (2020).

A broken cloud-filtering approach based on Gielen et al. (2014) is applied to the MAX-DOAS measurements to exclude MAX-DOAS aerosol and NO₂ scans influenced by the presence of clouds, which are known to potentially degrade the quality of the retrievals (Gielen et al., 2014; Wagner et al., 2014). Three sky conditions can be distinguished with this flagging approach: (1) clear sky, (2) homogeneous cloud coverage, and (3) broken cloud conditions. Retrievals under broken cloud conditions are rejected from the present study.

The profile retrieval was performed to estimate the mixing layer height of NO₂ (MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$). The MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is estimated per measurement scan, and it is the ratio of VCD${}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{main}}$ to the NO₂ near-surface concentration ($c_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{main}}$) as retrieved in the main azimuthal direction by the MMF inversion algorithm:

$$\begin{array}{}\text{(1)}& {\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{{\text{VCD}}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{main}}}{c_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{main}}}}.\end{array}$$

Therefore, during one measurement scan, two assumptions were made: (1) the homogeneous distribution of NO₂ inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ and (2) homogeneous MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ around the measurement site and its use in all the azimuthal directions. The validity of the second assumption is tested in Sect. 5.2.3.

5.2 Dual-scan MAX-DOAS retrieval method

A complete MAX-DOAS measurement scan is composed of two different sub-scans, as described in Sect. 2.1. The aerosol and NO₂ vertical profiles are retrieved from the elevation scan in the main azimuthal direction. In the other azimuthal directions, measurements are performed only in a single low elevation angle (2^∘), and therefore the retrieval of aerosol and NO₂ vertical profiles is not possible. Using the fact that the lowest elevation angles have the highest sensitivity to trace gases located near the surface due to the long light path in this layer, a new dual-scan MAX-DOAS retrieval strategy was developed here. This new retrieval strategy is an extension of the work presented in Dimitropoulou et al. (2020) and aims to retrieve the near-surface NO₂ box-averaged volume mixing ratios (VMRs) and the NO₂ VCDs at six different wavelengths. In Dimitropoulou et al. (2020), the applied dual-scan NO₂ MAX-DOAS retrieval was itself an adaptation of the parameterization technique proposed by Sinreich et al. (2013). More precisely, in the presence of sufficient aerosols in the atmosphere (i.e., sufficient aerosols to constrain the light path in a near-surface layer and ensure that the near-surface NO₂ concentration can be approximated by a near-surface box profile), the measured NO₂ dSCDs at one low elevation angle (2^∘) can be related to the near-surface NO₂ box-averaged concentration as follows:

$$\begin{array}{}\text{(2)}& {\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}=c_{{\mathrm{NO}}_{\mathrm{2}}}{L}_{{\mathrm{NO}}_{\mathrm{2}}},\end{array}$$

where dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is the differential slant column density of NO₂ and $c_{{\mathrm{NO}}_{\mathrm{2}}}$ its mean concentration along the differential effective light path, $L_{{\mathrm{NO}}_{\mathrm{2}}}$.

Consequently, the knowledge of the differential effective light-path length (i.e., $L_{{\mathrm{NO}}_{\mathrm{2}}}$) is crucial to derive the near-surface NO₂ concentrations. The oxygen collisional complex (O₄) can be used as a tracer for the effective light path in the atmosphere, as its concentration is well-known (it is the square of O₂ concentration). As a result, observed changes in the O₄ dSCDs can be directly attributed to changes in the light path due to the presence of particles like aerosols and clouds. $L_{{\mathrm{O}}_{\mathrm{4}}}$ is calculated by using Eq. (2) for O₄.

However, the direct use of the O₄ light-path length in the NO₂ retrieval is not possible under moderate to high pollution conditions, such as those in Brussels, because the profile shapes of O₄ and NO₂ are not the same. In Dimitropoulou et al. (2020), we used radiative transfer model (RTM) simulations to estimate a unitless correction factor, which accounts for these profile shape differences. This unitless correction factor indicates that under moderate to high pollution conditions, $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is equal to or smaller than $L_{{\mathrm{O}}_{\mathrm{4}}}$. For a correction factor equal or close to 1, $L_{{\mathrm{O}}_{\mathrm{4}}}$ is equal to $L_{{\mathrm{NO}}_{\mathrm{2}}}$, which means that there is a moderate to high aerosol load in the atmosphere during the measurement. On the other hand, correction factors smaller than unity are obtained for measurements performed under aerosol-free conditions or a thin MLH. Assuming a homogeneous NO₂ distribution inside the MLH, the MLH is derived from the NO₂ vertical profiles in the main azimuthal direction and is defined as the ratio of the NO₂ VCD to the near-surface concentration of NO₂. In Dimitropoulou et al. (2020), the RTM simulations were performed for eight different MLH values of aerosols and NO₂ in the range of 500–2000 m (i.e., eight different combinations) as well as for different measurement viewing geometries (solar zenith angle – SZA, relative azimuth angle – RAA, and the corresponding elevation angle of 2^∘). For every MAX-DOAS measurement, one value of the correction factor is given according to its viewing geometry and MLH value during the measurement. For further information, we refer the reader to Dimitropoulou et al. (2020).

In the present study, a new dual-scan NO₂ MAX-DOAS retrieval method, which is more suitable for interpreting multi-wavelength measurements than the previous approach (Dimitropoulou et al., 2020), is developed. It is presented in detail in the following subsection.

5.2.1 Developed dual-scan MAX-DOAS retrieval method

The main advantages of the new dual-scan NO₂ MAX-DOAS retrieval method (which are also the main differences with respect to Dimitropoulou et al., 2020) are the following: (1) the direct use of the measured O₄ dSCDs to estimate $L_{{\mathrm{NO}}_{\mathrm{2}}}$ for every measurement, (2) retrieval of near-surface aerosol extinction close to the ground, and (3) the exploitation of the wavelength dependency of the horizontal path representative of MAX-DOAS measurements for the retrieval of the horizontal distribution of aerosols (and therefore NO₂) around the measurement site. The latter is done using O₄ and NO₂ dSCDs measured at six different wavelengths. This new method is described below.

Assuming that the NO₂ vertical distribution can be approximated by a box profile of height equal to mixing layer height (MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$), the following equation can be used.

$$\begin{array}{}\text{(3)}& c_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{{\text{VCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}}}={\displaystyle \frac{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}\end{array}$$

This means that the NO₂ near-surface concentration can be expressed as a ratio of the dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ to the $L_{{\mathrm{NO}}_{\mathrm{2}}}$ (see Eq. 2) or as a ratio of the VCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ to the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$. When knowing MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ and VCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$, $L_{{\mathrm{NO}}_{\mathrm{2}}}$ can be simulated as follows.

$$\begin{array}{}\text{(4)}& L_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{simulated}}={\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{simulated}}{\displaystyle \frac{{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}}{{\text{VCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}}\end{array}$$

Here, NO₂ dSCDs and consequently $L_{{\mathrm{NO}}_{\mathrm{2}}}$ are simulated using the radiative transfer model VLIDORT version 2.7 (Spurr, 2006). Seasonal median MAX-DOAS NO₂ vertical profiles, as retrieved by applying the MMF inversion algorithm in the main azimuthal direction (see Sect. 5.1), show that the bulk (70 %) of the NO₂ concentration is located inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, which is expected since MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is estimated as the ratio of VCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ to the near-surface NO₂ concentration. On the other hand, this is not the case for aerosols (only 30 % of the aerosol content is seen to be located inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$). Considering this feature, for the VLIDORT simulations, the NO₂ a priori profiles are modeled as box profiles with a constant concentration equal to 1.5×10¹¹ molec. cm⁻³ from the surface to the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$. Two layers compose the aerosol a priori profiles: (1) the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ and (2) the free troposphere. The equation, which is applied to estimate the aerosol extinction profile a(z), is the following (see Wang et al., 2014):

$$\begin{array}{}\text{(5a)}& a\left(z\right)=\text{AOD}{\displaystyle \frac{p}{{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}}}\text{for}z\le {\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}\end{array}$$

and

$$\begin{array}{}\text{(5b)}& a\left(z\right)=b\left(\mathit{\xi },{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}},p\right)\exp \left(-{\displaystyle \frac{z}{\mathit{\xi }}}\right)\text{for}z>{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}},\end{array}$$

where AOD is the aerosol optical depth, p is the fraction of AOD inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, b is a normalizing constant for the exponential component (Wang et al., 2014), z is the simulation altitude grid, and ξ is the scaling height for the aerosols located outside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, which is set to 5 km (Wang et al., 2014). In the present study, the fraction of AOD located within the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is set to p=0.3 (see above). The effect of the p value and the NO₂ profile shape on the retrieved NO₂ near-surface VMRs and VCDs was investigated and considered in the error budget (see Sect. 5.2.2).

The MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is a known parameter, and it is estimated per measurement scan as the ratio of VCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ to the NO₂ near-surface concentration as retrieved in the main azimuthal direction by the MMF inversion algorithm.

The RTM simulations have nine input parameters in total, which are the elevation angle, SZA, RAA, AOD, MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, $c_{{\mathrm{NO}}_{\mathrm{2}}}$, AOD (p and ξ), and wavelength. It should be noted that the elevation angle is kept constant (i.e., 2^∘). For the six different wavelengths (343, 360, 380, 447, 477, and 530 nm), we separately perform RTM simulations, and $L_{{\mathrm{NO}}_{\mathrm{2}}}$ (see Eq. 4) is simulated for the assumed SZA, RAA, MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, $c_{{\mathrm{NO}}_{\mathrm{2}}}$, and AOD input scenarios presented in Table 2.

The simulated O₄ dSCDs are a function of the input parameter AOD. The relation between the simulated O₄ dSCDs and the input AOD values is shown in Fig. 4a. A piecewise cubic hermite interpolating polynomial fitting through the AOD as a function of the simulated O₄ dSCDs for each SZA, RAA, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ combination can be used in order to perform an inverse method (i.e., to estimate the near-surface aerosol extinction from the measured O₄ dSCDs). Additionally, in Fig. 4b, we can see that the relation between the O₄ dSCDs and the AOD values is valid for MAX-DOAS measurements.

Figure 4 (a) Dots: simulated AOD for NO₂ box profile of 1 km at 477 nm for an SZA of 40^∘ and RAA of 30^∘ as a function of the simulated O₄ dSCDs for the different AOD values (1, 0.8, 0.6, 0.4, 0.3, 0.1, and 0; see Table 2). Blue line: simulated AOD by applying an exponential fit through the data points. (b) Dots: estimated AOD in the six wavelengths used in the retrieval for 1 example day (11 September 2018). Red line: fitted line through the data points.

For every combination of all eight parameters (i.e., all the parameters of Table 2, except the AOD values), a polynomial fit of simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ as a function of simulated O₄ dSCDs is applied. Figure 5a shows simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ as a function of simulated O₄ dSCDs, and a second-order polynomial is fitted through the data points. We also observe in Fig. 5b, in which an example day of MAX-DOAS measurements is presented, that the $L_{{\mathrm{NO}}_{\mathrm{2}}}$ as a function of the measured O₄ dSCDs has the same relation as the simulated quantities. Since NO₂ is an optically thin absorber, $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is not a function of c${}_{{\mathrm{NO}}_{\mathrm{2}}}$, and consequently an $L_{{\mathrm{NO}}_{\mathrm{2}}}$ value can be estimated by using the measured O₄ dSCD for each measurement. Based on the corresponding SZA, RAA, measured O₄ dSCD, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, an $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is attributed to each low-elevation MAX-DOAS measurement through this polynomial fit. To express $L_{{\mathrm{NO}}_{\mathrm{2}}}$ as a function of four different parameters (i.e., O₄ dSCD, SZA, RAA, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$), $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is interpolated linearly at the O₄ dSCD, SZA, RAA, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ of each measurement. For example, a MAX-DOAS measurement with SZA =30^∘, RAA =60^∘, MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}=\mathrm{1}$ km, and measured O₄ dSCD $=\mathrm{6}\times {\mathrm{10}}^{\mathrm{43}}$ molec.² cm⁻⁵ will have an $L_{{\mathrm{NO}}_{\mathrm{2}}}$ equal to 15 km at 477 nm (see Fig. 5a).

Figure 5 (a) Dots (stars): simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ for NO₂ box profile of 1  km at 477 nm (530 nm) for an SZA of 30^∘ and RAA of 60^∘ as a function of the simulated O₄ dSCDs for the different AOD values (1, 0.8, 0.6, 0.4, 0.3, 0.1, and 0; see Table 2). Blue (red) line: second-order polynomial fit through the data points. (b) Dots: estimated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ in the six wavelengths used in the retrieval for 1 example day (11 September 2018). Red line: fitted line through the data points.

Based on this approach, the near-surface NO₂ concentration can be calculated at the six different wavelengths by using the measured dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ together with the simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ value (Eq. 2). The corresponding near-surface NO₂ VMRs are obtained by dividing the NO₂ concentrations by the air number density. To derive the air number density, we use monthly averaged pressure and temperature profiles over a 20-year period. These profiles are extracted from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis. In the last step, the tropospheric NO₂ VCD is calculated from the product of the near-surface NO₂ concentration with the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$.

Regarding aerosols, the AOD is estimated for every off-axis measurement (see Fig. 4a). The near-surface aerosol extinction is then calculated as the ratio between the aerosols inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ (i.e., AOD times p) and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$. The near-surface aerosol extinction refers to the layer that extends from the surface to the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$. As mentioned above, around 30 % of the total aerosols are expected to be found inside this layer.

The effect of SZA, RAA, and MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ on the simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is investigated in the Supplement. First, the simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is presented in Fig. S1 as a function of RAA for different MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ and wavelengths as well as a single AOD and SZA value. $L_{{\mathrm{NO}}_{\mathrm{2}}}$ strongly depends on MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$. The lower the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, the shorter the $L_{{\mathrm{NO}}_{\mathrm{2}}}$. The same NO₂ concentration and aerosol load are used for the three different MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ scenarios. So, when aerosols are concentrated in a thin layer (i.e., MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}=\mathrm{0.5}$ km), $L_{{\mathrm{NO}}_{\mathrm{2}}}$ becomes shorter. Secondly, we observe that $L_{{\mathrm{NO}}_{\mathrm{2}}}$ depends on RAA. The larger the RAA, the longer the $L_{{\mathrm{NO}}_{\mathrm{2}}}$. In Fig. S2, simulated $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is plotted for each wavelength and each considered MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ as a function of SZA (at a constant AOD and RAA). $L_{{\mathrm{NO}}_{\mathrm{2}}}$ depends strongly on SZA. The highest dependency is observed for large SZA values at which $L_{{\mathrm{NO}}_{\mathrm{2}}}$ becomes maximum. Finally, in both Figs. S1 and S2, we observe that $L_{{\mathrm{NO}}_{\mathrm{2}}}$ becomes longer with wavelength, which is expected because of the less pronounced Rayleigh and Mie scattering at longer wavelengths.

An example of dual-scan MAX-DOAS retrieval is shown in Fig. 6. Based on the RTM simulations described above, $L_{{\mathrm{NO}}_{\mathrm{2}}}$ is derived for the wavelengths of interest, and ultimately, near-surface NO₂ concentrations and tropospheric NO₂ VCDs are estimated. In the last step, the near-surface aerosol extinction values are assigned to the six different wavelengths.

Figure 6 (a) Corresponding $L_{{\mathrm{NO}}_{\mathrm{2}}}$, (b) near-surface NO₂ concentrations, (c) NO₂ VCDs, and (d) aerosol optical densities as a function of the six wavelengths used in the retrieval (11 September 2018, 11:51 UTC, 123.5^∘ azimuthal direction).

5.2.2 Uncertainty budget

To estimate uncertainties in the dual-scan parameterized NO₂ near-surface concentration and VCD, the standard error propagation method is used as

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{c_{{\mathrm{NO}}_{\mathrm{2}}}}^{\mathrm{2}}={\left({\mathit{\sigma }}_{{\mathrm{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{\displaystyle \frac{\partial c_{{\mathrm{NO}}_{\mathrm{2}}}}{\partial {\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{L_{{\mathrm{NO}}_{\mathrm{2}}}}{\displaystyle \frac{\partial c_{{\mathrm{NO}}_{\mathrm{2}}}}{\partial L_{{\mathrm{NO}}_{\mathrm{2}}}}}\right)}^{\mathrm{2}},\end{array}$$

which is solved as

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{c_{{\mathrm{NO}}_{\mathrm{2}}}}^{\mathrm{2}}={\left({\displaystyle \frac{{\mathit{\sigma }}_{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}\right)}^{\mathrm{2}}+{\left(-{\displaystyle \frac{{\mathit{\sigma }}_{L_{{\mathrm{NO}}_{\mathrm{2}}}}{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{L_{{\mathrm{NO}}_{\mathrm{2}}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}.\end{array}$$

Then,

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{c_{{\mathrm{NO}}_{\mathrm{2}}}}^{\mathrm{2}}={\left({\mathit{\sigma }}_{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{\displaystyle \frac{c_{{\mathrm{NO}}_{\mathrm{2}}}}{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{L_{{\mathrm{NO}}_{\mathrm{2}}}}{\displaystyle \frac{c_{{\mathrm{NO}}_{\mathrm{2}}}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}\right)}^{\mathrm{2}},\end{array}$$

where

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{\sigma }}_{L_{{\mathrm{NO}}_{\mathrm{2}}}}^{\mathrm{2}}=& {\left({\displaystyle \frac{L_{{\mathrm{NO}}_{\mathrm{2}}}}{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}\left(\text{sim}\right)}}{\mathit{\sigma }}_{{\text{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{L_{{\mathrm{NO}}_{\mathrm{2}}}}{{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}}}{\mathit{\sigma }}_{{\text{MLH}}_{{\mathrm{NO}}_{\mathrm{2}}}}\right)}^{\mathrm{2}}.\end{array}\end{array}$$

According to Kreher et al. (2020) and Bösch et al. (2018), in urban or suburban polluted conditions, the use of the DOAS fit uncertainty of NO₂ for the dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ uncertainty is not appropriate because the dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ uncertainty is mostly driven by atmospheric variability as well as spatial and temporal fluctuations in the O₄ and NO₂ fields. In this study, a conservative value of 3.5×10¹⁵ molec. cm⁻² is attributed to ${\mathit{\sigma }}_{{\mathrm{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}$ (Kreher et al., 2020). This represents an error of up to 6.0 % in the NO₂ dSCDs in the visible range (477 nm).

The second error source is related to the estimation of $L_{{\mathrm{NO}}_{\mathrm{2}}}$ from the RTM simulations. To estimate this error, sensitivity tests on the input aerosol and NO₂ vertical profiles were performed. The fraction of aerosols located inside the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ (40 % and 60 % instead of 30 %) and the NO₂ profile shape (linearly decreasing instead of box) were modified. The error related to the RTM simulations of the simulated dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is about 9.6 % in the visible range (477 nm). Additionally, according to Dimitropoulou et al. (2020), the uncertainty related to MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ is about 4 % in the visible range (477 nm).

Combining all the error sources, the total uncertainties in the NO₂ near-surface concentration are about 12.7 %, 12.4 %, 11.4 %, 11.1 %, 11.3 %, and 12.6 % at 343, 360, 380, 447, 477, and 530 nm, respectively.

Finally, the total uncertainties in the NO₂ VCD are about 12.0 %, 11.7 %, 10.7 %, 10.3 %, 10.6 %, and 11.9 % at 343, 360, 380, 447, 477, and 530 nm, respectively.

5.2.3 Verification of the dual-scan MAX-DOAS retrieval method

The consistency check and verification of the dual-scan MAX-DOAS retrieval method in Uccle are based on two different correlative comparisons.

The consistency check compares the NO₂ near-surface VMRs and tropospheric VCDs retrieved by the dual-scan parameterization in the main azimuthal direction to the same quantities retrieved with the MMF inversion algorithm at the two main wavelengths (360 and 477 nm). As can be seen in Figs. 7 and 8, both datasets are in good agreement, with correlation coefficient values in the range of 0.86 to 0.95 and slope values close to unity for all four comparisons.

Figure 7 Comparison between MMF and parameterized NO₂ near-surface VMR at 477 nm (visible, a) and 360 nm (UV, b), as derived from the main azimuthal direction (i.e., 35.5^∘ azimuthal direction).

Figure 8 Comparison between MMF and parameterized NO₂ VCD at 477 nm (visible, a) and 360 nm (UV, b), as derived from the main azimuthal direction (i.e., 35.5^∘ azimuthal direction).

The verification step is based on the same type of comparison as the first one but for three additional azimuthal directions in which elevation scans, and hence profile retrievals, are available for some periods. Onward from 3 July 2019, elevation scans were performed in these three additional azimuthal directions to complement the already existing measurement setup. These elevation scans were performed once per day around noon in the 11, 105, and 262.5^∘ azimuthal directions. Figure 9 shows the comparison between near-surface NO₂ VMRs and tropospheric VCDs retrieved by the dual-scan parameterization method and the corresponding results obtained with the MMF inversion algorithm. Overall good agreement is obtained (R=0.79 and 0.84 for near-surface VMR and VCD, respectively). We observe that the comparison concerning the near-surface NO₂ VMR seems to be noisier than in the main azimuth direction. This is mainly due to the use of the MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ calculated in the main azimuthal direction for all the different azimuth angles in the dual-scan method. Additionally, the parameterization technique slightly underestimated the near-surface NO₂ VMR (s=0.84), while a slope value of 1.00 is obtained for tropospheric VCDs.

Figure 9 Visible range: comparison between MMF and (a) parameterized NO₂ near-surface VMR as well as (b) parameterized NO₂ VCD at three different azimuthal directions, as indicated in the color bar (11, 105, and 262.5^∘ azimuthal directions). The elevation scans in these azimuthal directions were performed once per day from 3 July 2019.

5.3 Horizontal distribution inversion approach

The parameterized NO₂ near-surface concentrations at the six different wavelengths are used as input in a new horizontal distribution inversion approach. As parameterized NO₂ near-surface concentrations, we refer to the conversion of the measured NO₂ dSCDs (i.e., at the elevation angle of 2^∘) to near-surface NO₂ concentrations by applying the dual-scan MAX-DOAS retrieval method as described in Sect. 5.2. Figure 10 shows a sketch of the assumed horizontal box model configuration, in which successive boxes of concentration c_N between the horizontal distances d_N−1 and d_N from the MAX-DOAS instrument are considered along the light path. The index N is equal to the total number of successive boxes.

Figure 10 Schematic representation of the six different $L_{{\mathrm{NO}}_{\mathrm{2}}}$ (i.e., one horizontal line for each wavelength) used in the new horizontal distribution inversion approach. The length of each line shows the sensitivity of each wavelength as a function of the horizontal distance. The shortest line represents the smallest wavelength.

The different horizontal lines illustrate the horizontal extent (or differential effective light path as described in Sect. 5.2) in which the NO₂ near-surface concentrations are extended for the six different wavelengths. Generally, the MAX-DOAS horizontal sensitivities are longer for larger wavelengths because of the less pronounced Rayleigh scattering (see also Fig. 6; Ortega et al., 2016; Dimitropoulou et al., 2020). In Fig. 10, the shortest line represents the smallest wavelength's horizontal sensitivity (343 nm) and the longest line the largest wavelength's horizontal sensitivity (530 nm). As can be seen in the sketch, the effective horizontal light path at the six different wavelengths passes through a different number of horizontal bins.

The parameterized NO₂ near-surface concentrations at the different wavelengths are the mean concentrations along the horizontal effective light paths (see Sect. 5.2), which are also called differential effective light paths because they are linked to the dSCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$. When having information coming from one wavelength only, it is not possible to know how the NO₂ is distributed along this light path. In the present work, the knowledge of mean NO₂ concentrations at six different wavelengths is used to retrieve a horizontal NO₂ profile, assuming the horizontal box model described in Fig. 10. This new retrieval method is described below.

The measurement vector y consists of the six retrieved surface concentrations (called ${\stackrel{\mathrm{‾}}{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$; see the method presented in Sect. 5.2) at the six different wavelengths. These near-surface concentrations can be expressed as functions of the different effective light paths ($L_{{\mathrm{NO}}_{\mathrm{2}}}$) and correspond to the average surface concentrations along those $L_{{\mathrm{NO}}_{\mathrm{2}}}$.

$$\begin{array}{}\text{(10)}& \mathit{y}={\mathit{F}}_{\mathrm{meas}}\left({\mathit{c}}_{{{\mathrm{NO}}_{\mathrm{2}}}_{\mathrm{true}}}\right)={\stackrel{\mathrm{‾}}{\mathit{c}}}_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{{\mathrm{dSCD}}_{{\mathrm{NO}}_{\mathrm{2}}}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}\end{array}$$

F_calcul, which represents the forward model, can be expressed as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{F}}_{\mathrm{calcul}}\left({\mathit{c}}_{{{\mathrm{NO}}_{\mathrm{2}}}_{\mathrm{true}}}\right)=& {{\displaystyle \frac{\mathrm{1}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}}^{\mathrm{2}}{\int }_{\mathrm{0}}^{L_{{\mathrm{NO}}_{\mathrm{2}}}}{\mathit{e}}^{-\frac{\mathit{x}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}{\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}\left(x\right)\\ & \cdot {\displaystyle \frac{\mathrm{d}x}{\cos \left(a\right)}},\end{array}\end{array}$$

where x is the horizontal distance, a is the elevation angle of the measurement, and ${\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$ is NO₂ near-surface concentration as a function of x, the distance from the MAX-DOAS instrument.

Our retrieval of the horizontal distribution of ${\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$ is based on the inversion theory (Rodgers, 2000), in which a horizontal profile ${\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$ (state vector) is retrieved given an a priori horizontal profile x_α, the measurement vector y, the matrix of the weighting function K, the uncertainty covariance matrix of the a priori S_α and the uncertainty covariance matrix of the measurement S_e.

$$\begin{array}{}\text{(12)}& {\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}={\mathit{x}}_{\mathrm{a}}+{\left({\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}(\mathit{y}-\mathbf{K}{\mathit{x}}_{\mathrm{a}})\end{array}$$

The weighting function indicates the sensitivity of the measurement vector to a change in the horizontal profile. Following the Beer–Lambert law, the weighting functions follow an exponential decrease between the instrument and $L_{{\mathrm{NO}}_{\mathrm{2}}}$. Additionally, we should consider that the light path through the horizontal boxes is slanted and not horizontal by including the cosine of the elevation angle (i.e., 2^∘ elevation angle). The weighting functions are given by the following equation.

$$\begin{array}{}\text{(13)}& \mathbf{K}\left(\mathit{\lambda },x\right)={\displaystyle \frac{\mathrm{1}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}{\mathit{e}}^{-\frac{\mathit{x}}{L_{{\mathrm{NO}}_{\mathrm{2}}}}}{\displaystyle \frac{\mathrm{\Delta }x}{\cos \left(a\right)}}\end{array}$$

An example of weighting functions is presented in Fig. 11. The sensitivity decreases exponentially up to a distance corresponding to the differential effective light-path length of each measurement. More precisely, each measurement is highly sensitive to the MAX-DOAS instrument location. This sensitivity decreases exponentially as a function of the horizontal distance. Then, it reaches a value equal to $\mathrm{1}/e$ to the horizontal distance equal to the differential effective light-path length of each measurement. It should be noted that since NO₂ is an optically thin absorber, the measurements depend linearly on each horizontal box's concentration. For this reason, OEM for the linear case is considered here, and only one inversion step is needed (see Eq. 12).

Figure 11 Examples of weighting functions used in the new horizontal distribution inversion approach (11 September 2018, 11:26 UTC).

The selected output horizontal grid for the retrieval extends from the MAX-DOAS instrument to the horizontal distance in which the weighting function of the largest wavelength (i.e., 530 nm) decreases at 10 % of its value at the MAX-DOAS instrument's location per azimuthal direction and consists of successive boxes of 0.5 km thickness on the horizontal axis.

Since this inversion problem is ill-conditioned, more than one horizontal NO₂ profile can be consistent with the measurement vector. To reject unrealistic solutions, the a priori profile x_a and its uncertainty covariance matrix must be included in the retrieval. In the OEM, the a priori information usually comes from an independent source, like a model or other correlative measurements. In the present study, RIO model data were chosen as a priori. Seasonal average maps of RIO NO₂ near-surface concentration are constructed (see Fig. S3), and then seasonal averages of RIO NO₂ near-surface concentration horizontal profiles were calculated in each azimuthal direction and interpolated on the retrieval's horizontal grid by regridding the initial 4×4 km² spatial resolution to a finer one (see Fig. 12). The shape of the RIO a priori NO₂ profiles per azimuthal direction stays the same during different seasons of the year, indicating that the wind effect on NO₂ transportation disappears by the seasonal averaging and that the same sources contribute to the NO₂ horizontal field. A mean scaling factor equal to the mean ratio between the measured and RIO NO₂ near-surface concentrations is applied because of the systematic underestimation of NO₂ near-surface concentrations by MAX-DOAS when compared to in situ measurements (see Dimitropoulou et al., 2020, and Sect. S1).

Figure 12 Example of seasonal RIO a priori NO₂ horizontal profiles for the new horizontal distribution inversion approach as a function of the horizontal distance from the MAX-DOAS instrument in six different azimuthal viewing directions, before the application of the scaling factor.

For the aerosol horizontal distribution retrieval, there are no sufficient independent measurements that provide information about the horizontal distribution of AOD and can serve as an a priori AOD profile. Therefore, a constant a priori AOD profile is used in the AOD retrieval based on CIMEL observations. An AOD equal to 0.18, which is the yearly averaged AOD value from CIMEL at 477 nm, is used. To construct the near-surface aerosol extinction a priori profiles, 30 % of the total amount of AOD is considered to be located inside the MLH (i.e., known for each MAX-DOAS vertical scan from the MMF inversion algorithm; see Sect. 5.1).

The diagonal elements of the S_a matrix are set equal to the square of a scaling factor times the NO₂ concentration a priori profile. The non-diagonal elements, which account for correlation between the different horizontal grid cells, are set as follows (Barret et al., 2002):

$$\begin{array}{}\text{(14)}& {\mathbf{S}}_{{\mathrm{a}}_{ij}}=\sqrt{{\mathbf{S}}_{{\mathrm{a}}_{ii}}{\mathbf{S}}_{{\mathrm{a}}_{ij}}\exp \left(-\ln \left(\mathrm{2}\right){\left({\displaystyle \frac{x_i-x_j}{\mathit{\gamma }}}\right)}^{\mathrm{2}}\right)},\end{array}$$

where x_i and x_j are the horizontal distances at the ith and jth horizontal boxes, and γ is half of the correlation length. For NO₂, γ is set equal to 3.5 km, and for aerosols, γ is set equal to 2 km. To eliminate inversion instabilities, S_a elements which are smaller than 0.1 % of the maximum S_a element are set equal to zero.

To estimate the correlation length, a covariance matrix was constructed by exploiting the airborne observations above Brussels (28 June 2019). The airborne observations have a spatial resolution of approximately 100×100 m². NO₂ horizontal profiles were constructed in different azimuthal directions at a spatial resolution of 500×500 m², expanding from the MAX-DOAS position to a maximum distance of 20 km, and were used to calculate a covariance matrix. A correlation length equal to 7 km and, consequently, a gamma value equal to 3.5 km are found to be representative for the NO₂ horizontal profiles in Brussels.

Using this correlation length, a variance of 45 % is used. This choice was made based on the seasonal variance of the RIO a priori profiles compared to their seasonal mean value. It is found that the seasonal variance of RIO observations has a mean value of 45 %. Additionally, it is found to be a good compromise for obtaining reasonable retrieval results, e.g., in terms of information content, while avoiding unrealistic oscillations in the retrieved aerosol and NO₂ profiles. It should be noted that for aerosols, a variance of 65 % is used.

The measurement covariance matrix S_e is chosen to be diagonal, with elements corresponding to the uncertainties of the dual-scan parameterized NO₂ near-surface concentration (see Sect. 5.2.2).

Three examples of the retrieved NO₂ horizontal profile are presented in Fig. 13, together with corresponding measured and simulated ${\stackrel{\mathrm{‾}}{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$ at the six different wavelengths for 2 July 2018 (Fig. 13a; low NO₂ abundance condition), 11 September 2018 (Fig. 13b; medium NO₂ abundance condition), and 30 September 2018 (Fig. 13c; medium NO₂ abundance condition). The root mean square (rms) is calculated between measured and simulated NO₂ near-surface concentrations of the horizontal retrieval normalized by the mean of the measured NO₂ near-surface concentrations (upper panels in Fig. 13). As the NO₂ values become larger, the agreement between measured and simulated ${\stackrel{\mathrm{‾}}{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$, expressed via the rms value, is improved.

Figure 13 Measured and retrieved NO₂ near-surface concentrations (upper panels in a, b, and c) at the six different wavelengths (i.e., horizontal distances) as a function of the estimated horizontal distances and (lower panels in a, b, and c) the retrieved NO₂ near-surface horizontal profile and a priori profile. The panels correspond to (a) 11 September 2018 at 12:91 UTC with a 265^∘ azimuthal direction, (b) 2 July 2018 at 10:42 UTC with a 25^∘ azimuthal direction, and (c) 30 September 2018 at 07:96 UTC with a 167.5^∘ azimuthal direction.

Similarly, examples of measured and retrieved near-surface aerosol extinction coefficients and retrieved aerosol horizontal profiles are shown in Fig. 14 for different aerosol load conditions (low, Fig. 14a; medium, Fig. 14b; and high, Fig. 14c) over the Brussels-Capital Region. We observe that the agreement between simulated and measured near-surface aerosol extinction coefficients at the six different wavelength tends to be worse than for NO₂. This could be due to the use of a constructed (constant) a priori aerosol horizontal profile due to the lack of information on the aerosol extinction horizontal distribution in the Brussels-Capital region.

Figure 14 Measured and retrieved near-surface aerosol extinction (upper panels in a, b, and c) at the six different wavelengths (i.e., horizontal distances) as a function of the estimated horizontal distances and (lower panels in a, b, and c) the retrieved near-surface aerosol extinction horizontal profile. The panels correspond to (a) 11 September 2018 at 11:48 UTC with a 167.5^∘ azimuthal direction, (b) 20 September 2018 at 07:08 UTC with a 25^∘ azimuthal direction, and (c) 2 June 2018 at 18:16 UTC with a 35.5^∘ azimuthal direction.

An essential condition of the dual-scan MAX-DOAS retrieval and the new horizontal inversion approach at six different wavelengths is the increasing trend of the horizontal sensitivity as a function of wavelength. Consequently, every wavelength is sensitive to a different horizontal region, and the six different wavelengths can be used to retrieve the horizontal distribution of aerosols and trace gases. Sensitivity tests were conducted in which simulated $L_{{\mathrm{O}}_{\mathrm{4}}}$ is expressed as a function of the six different wavelengths for different aerosol conditions. As can be seen in Fig. S4, the linear relationship between $L_{{\mathrm{O}}_{\mathrm{4}}}$ (and $L_{{\mathrm{NO}}_{\mathrm{2}}}$) and wavelength exists for AOD values ranging from 0 to 1. An AOD equal to unity is chosen as the maximum AOD of the simulations because in Uccle AOD values rarely exceed 1 (see in https://aeronet.gsfc.nasa.gov/, last access: 11 July 2022, for the Brussels measurement site). Therefore, the relation stays linear as the aerosol load changes for the conditions observed in Uccle. The only condition leading to nonlinearity is when clouds are present. However, as explained in Sect. 5.1, a cloud-filtering approach is applied, rejecting the broken cloud scenes, which are the more problematic ones.

5.4 Characterization of the retrieval

To characterize the retrieval, the averaging kernels, AK, play a crucial role. The AK matrix is calculated as follows (Rodgers, 2000).

$$\begin{array}{}\text{(15)}& \mathbf{AK}={\displaystyle \frac{\mathrm{d}{\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}}{\mathrm{d}{\mathit{c}}_{{{\mathrm{NO}}_{\mathrm{2}}}_{\mathrm{true}}}}}={\left({\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathbf{K}\end{array}$$

The AKs are the rows of the AK matrix. They present the sensitivity of the retrieved (${\mathit{c}}_{{\mathrm{NO}}_{\mathrm{2}}}$) to the true (${\mathit{c}}_{{{\mathrm{NO}}_{\mathrm{2}}}_{\mathrm{true}}}$) atmospheric profile. Ideally, the AK matrix should be an identity matrix. In Fig. 15, an example of selected AKs is shown. As can be seen that, for distances smaller than the first measurement (e.g., near-surface NO₂ concentration retrieved at 343 nm), the AKs have maximum values and decrease exponentially as a function of the horizontal distance from the instrument. For this particular example, the AKs do not exceed the values of 0.12.

Figure 15 Example of NO₂ averaging kernels. They are calculated for observations on 11 September 2018 at 11:51 UTC and 300^∘ azimuthal direction.

Another important point about the retrieval is the trace of the AK matrix, which refers to the number of degrees of freedom for signal (DOFS). The DOFS are an indication of the number of independent pieces of information that one can retrieve from the measurements. Ideally, the DOFS would be equal to the number of horizontal boxes for the horizontal distribution. In reality, the DOFS are lower because of the limited horizontal resolution of the measurements. In Fig. 15, the DOFS are close to 2, which means that two independent pieces of information are contained in the measurements for this particular example.

In the present work, the total retrieval error is equal to the error related to the measurement noise and the smoothing of the true atmospheric profile. According to Rodgers (2000), the noise covariance matrix is estimated as

$$\begin{array}{}\text{(16)}& {\mathbf{S}}_{\mathrm{meas}}={\mathbf{GS}}_{\mathrm{e}}{\mathbf{G}}^{\mathrm{T}},\end{array}$$

with G being the gain matrix:

$$\begin{array}{}\text{(17)}& \mathbf{G}={\left({\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}.\end{array}$$

Then, the retrieval noise error is given as the square root of the diagonal elements of the noise covariance matrix.

The smoothing error is calculated as follows:

$$\begin{array}{}\text{(18)}& {\mathbf{S}}_{\mathrm{smooth}}=\left(\mathbf{AK}-\mathbf{I}\right){\mathbf{S}}_x{\left(\mathbf{AK}-\mathbf{I}\right)}^{\mathrm{T}},\end{array}$$

with S_x being a realistic covariance metric of the true NO₂ horizontal profile.

The horizontal profiles of the measurement and smoothing error in percentage are shown in Fig. 16. As can be seen, the smoothing error is significantly larger than the measurement error (range of 3 %–10 % and 14 %–40 %, respectively). The smoothing error also becomes larger as the horizontal distances from the instrument become larger. This is mainly because of the exponential decrease in the sensitivity as a function of the horizontal distance (see weighting functions in Fig. 11) and, consequently, the larger impact of the difference between the a priori profile and the true state of the atmosphere.

Figure 16 Example of the NO₂ measurement and smoothing error in percentage for 2 July 2018 at 10:42 UTC and 25^∘ azimuthal direction.

To eliminate the unsuccessful retrievals, the percentage of accepted retrievals with respect to the total number of retrievals during the four seasons is investigated when a specific filtering on rms and DOFS is applied (see Table 3 and Figs. S7 and S8). As we can see in Fig. S8, DOFS are in the range of 1.0–1.6. From these tests, it is found that most of the retrievals have DOFS larger than 1.2 (see Fig. S8). The rms is defined as the root mean square deviation between measured and simulated $c_{{\mathrm{NO}}_{\mathrm{2}}}$ normalized by the mean of the measured $c_{{\mathrm{NO}}_{\mathrm{2}}}$ (e.g., same rms as in Fig. 13). Table 3 and Fig. S7 indicate that rms values are in the range of 0 %–30 %, and most of the retrievals have an rms smaller than 9 % with a median rms value of around 5.6 % during all seasons. Based on these investigations, DOFS >1.2 and rms <9 % are used as retrieval quality control criteria.

Table 3 Seasonally averaged root mean square (rms) and DOFS values. The rms is calculated between measured and retrieved NO₂ near-surface concentrations of the horizontal retrieval (Fig. 13). DOFS represent the degrees of freedom of the horizontal retrieval (Fig. 15). The percentage of the accepted retrievals is presented for the different selection criteria.

6 Retrieval results and discussion

6.1 Example of daily horizontal NO₂ distribution

The variation of the MAX-DOAS horizontal distribution of tropospheric NO₂ VCDs as a function of time over the course of 28 June 2019 is presented in Fig. 17. This particular day is chosen because airborne measurements took place above the Brussels region (see Sect. 6.2). The horizontal NO₂ profiles are plotted per azimuthal direction with the horizontal axis showing the time in UTC and the vertical axis the horizontal distance in kilometers. Because of the quality check on the retrieved NO₂ horizontal profiles (see Sect. 5.3), some profiles were rejected (e.g., azimuthal direction equal to 262.5 and 265^∘).

Figure 17 Diurnal variation of the retrieved NO₂ horizontal profiles per azimuthal direction as a function of time (UTC) for 28 June 2019.

During this day, maximum NO₂ columns are mainly observed around 05:00 and 10:00 UTC, which correspond to 07:00 and noon local time. Early in the morning (05:00 UTC), high NO₂ columns are expected to be observed because of the low MLH (MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$ in the range of 300–600 m height) in combination with the morning rush-hour NO₂ emissions. Around 10:00 UTC, the maximum NO₂ columns are detected in the north (N), northeast (NE), and northwest (NW) direction (see Fig. 18). The same NO₂ horizontal distribution is found when investigating the NO₂ near-surface concentrations for this day (see Fig. S5). In the Brussels region, the main emission sources are located in the N and west (W) parts of the city and are linked to the motorway around Brussels (the so-called Ring), the Brussels city center, and the Drogenbos power plant (NW direction). Concerning the NO₂ peaks, they are located at a distance of around 0 to 8 km from the measurement site. It can be seen from Fig. 18 that the Ring, the Brussels city center, and the Drogenbos power plant are located within these distances. As measured by the meteorological station on the BIRA-IASB rooftop, the wind was coming from the E direction during that day, resulting in the progressive displacement of the NO₂ peak from the NNE to the W direction. On the contrary, the azimuthal directions pointing towards a large forested area (i.e., 62.5, 75, and 105^∘), the Bois de la Cambre, detect considerably lower NO₂ columns than the other directions.

Figure 18 Maps of hourly averaged NO₂ horizontal profiles per azimuthal direction for 28 June 2019 corresponding to Fig. 17. The wind direction is shown with the black arrow. The black square shows the MAX-DOAS instrument location, the black polygon the national airport, the black dots the NO₂ hotspots emitting more than 10 kg NO_x h⁻¹ (Emission Inventory of the Belgian Interregional Environment Agency, 2017; Frederik Tack, personal communication, 2019), and the black line the Brussels Ring road.

Maximum near-surface aerosol extinction coefficient values are observed all day long and detected in the N, NW, and NE direction (see Fig. S6). NO₂ and aerosol peaks are co-located towards the N and NW direction.

6.2 MAX-DOAS horizontal NO₂ distribution versus airborne, car-mobile DOAS, and TROPOMI: 28 June 2019 case study

The APEX tropospheric NO₂ columns are compared to the tropospheric NO₂ horizontal distribution as retrieved by applying our new mapping MAX-DOAS technique to the 28 June 2019 measurements. During the same day, TROPOMI pixels (OFFL 010302 product; see Table 1) selected over the Brussels region are compared to MAX-DOAS observations. During this day, the TROPOMI overpass time was at 12:19 UTC. MAX-DOAS horizontal profiles of tropospheric NO₂ VCDs are selected around the TROPOMI overpass time (±1 h). The horizontal profile of MAX-DOAS NO₂ VCDs for each horizontal line of sight has a horizontal sampling of 0.5 km (see Fig. 13). The MAX-DOAS NO₂ VCDs for the horizontal segment crossing a TROPOMI pixel and located inside the pixel are averaged and compared to the corresponding TROPOMI NO₂ VCD. It should be noted that the MAX-DOAS segments are weighted by their relative length inside each pixel. APEX observations located inside each TROPOMI pixel were used to assign one APEX NO₂ VCD value per pixel. Maps of co-located TROPOMI, AEROMOBIL, averaged MAX-DOAS, and averaged APEX NO₂ VCDs for 28 June 2019 are shown in Fig. 19. Two maps of APEX and MAX-DOAS observations are presented: one with APEX and MAX-DOAS in their initial resolution and one with spatially averaged APEX and MAX-DOAS observations in the area covered by a TROPOMI pixel. The NO₂ plume as detected by APEX is covering the NW, N, and NE parts of the Brussels region. MAX-DOAS successfully detected the same NO₂ plume in the NW and N but not in the NE direction. The correlation and agreement between APEX and MAX-DOAS observations are very good (R=0.89 and s=1.22).

Figure 19 (a) Tropospheric NO₂ VCD as detected by the APEX instrument in its initial spatial resolution. Tropospheric NO₂ VCD maps (TROPOMI pixels) as retrieved over Brussels on 28 June 2019 by (b) APEX, (d) MAX-DOAS, and (f) TROPOMI (overpass time at 12:19 UTC). Tropospheric NO₂ VCD as retrieved by (c) MAX-DOAS and (e) AEROMOBIL (between 08:30 and 15:42 UTC) in its initial spatial resolution. The black square shows the MAX-DOAS instrument location, the black polygon the national airport, the black dots the NO₂ hotspots emitting more than 10 kg NO_x h⁻¹ (Emission Inventory of the Belgian Interregional Environment Agency, 2017; Frederik Tack, personal communication, 2019), and the black line the Brussels Ring road.

During the S5PVAL-BE flight over Brussels, the AEROMOBIL was used to measure the spatial distribution of tropospheric NO₂ columns, mainly over the Ring road of Brussels, with a measurement start time at 08:30 UTC and end time at 15:42 UTC. The AEROMOBIL NO₂ VCDs, which are located closest in both time and space to each MAX-DOAS horizontal grid, are compared to the corresponding MAX-DOAS VCDs (see Fig. 19c and e). AEROMOBIL and MAX-DOAS agree perfectly on the location of maximum (i.e., NW direction) and minimum (i.e., SE direction) NO₂ tropospheric VCDs (Fig. 19c and e). We can observe in Fig. 20b that the correlation coefficient is moderate (R equal to 0.79) and the slope value is equal to 0.59. The correlation plot between the two datasets reveals that AEROMOBIL gives higher NO₂ tropospheric VCDs compared to MAX-DOAS ones. This finding could be partly explained by the fact that AEROMOBIL follows busy routes on which the NO₂ tropospheric VCDs reach maximum values because of the contribution of NO₂ production from vehicle engines via fossil fuel combustion.

During the TROPOMI overpass (i.e., 12:19 UTC) above the Brussels-Capital Region, dual-scan MAX-DOAS tropospheric NO₂ columns are retrieved, as can be seen in Fig. 18. The correlation between TROPOMI and MAX-DOAS tropospheric NO₂ columns during the day of the airborne measurements above Brussels is presented in Fig. 20c. Excellent agreement is obtained, with a correlation coefficient value equal to 0.88. The slope value is equal to 0.70. During that day, MAX-DOAS and TROPOMI are in good agreement but TROPOMI tends to underestimate the tropospheric NO₂ columns. It should be noted that during that day, the range of observed NO₂ VCDs is from 4.4×10¹⁵ to 9.5×10¹⁵ molec. cm⁻², as retrieved by the MAX-DOAS observations.

Figure 20 Scatter plot between (a) the tropospheric NO₂ columns derived by airborne measurements (APEX) and the MAX-DOAS observations, (b) the tropospheric NO₂ columns derived by car-mobile DOAS measurements (AEROMOBIL) and the MAX-DOAS observations, and (c) the tropospheric NO₂ columns derived by MAX-DOAS observations and the TROPOMI tropospheric NO₂ columns over Brussels on 28 June 2019.

6.3 Comparison between MAX-DOAS horizontal NO₂ distribution and TROPOMI observations

To compare the TROPOMI and MAX-DOAS tropospheric NO₂ columns, a similar approach is used as in Sect. 6.2. Additionally, TROPOMI and MAX-DOAS tropospheric NO₂ columns are compared on a seasonal basis, and seasonally averaged maps of those VCDs in the area covered by the TROPOMI pixels are created. To generate these maps, the ensemble of TROPOMI pixels recorded on 28 June 2019 is chosen as a reference and TROPOMI pixels that coincide with this reference grid are averaged. The daily horizontal profiles of MAX-DOAS NO₂ columns are averaged on the daily TROPOMI grids and then the reference grid is used to create the seasonally averaged MAX-DOAS maps.

The seasonally and annually averaged maps of TROPOMI and MAX-DOAS NO₂ VCDs are presented in Figs. 21 and 22. Only pixels including at least 20 comparison days are taken into account in the analysis. It is found that the locations of the NO₂ peaks and dips show a reasonably high degree of similarity between TROPOMI and MAX-DOAS during all seasons, except summer. The NO₂ peaks appear mainly above Brussels city center, the Drogenbos power plant (W direction), and the NW part of the Ring road, which are the main known emission sources, as mentioned earlier. These maps also indicate that the tropospheric NO₂ column over the Brussels area has a clear seasonal cycle, with a maximum during winter.

Figure 21 Seasonal tropospheric NO₂ VCD grids (TROPOMI grids) as retrieved over Brussels by TROPOMI and MAX-DOAS. The black square shows the MAX-DOAS position, the black polygon the national airport, the black dots the NO₂ hotspots, and the black line the Brussels Ring motorway.

Figure 22 Annual (e.g., based over the 2 years of observations) tropospheric NO₂ VCD grids (TROPOMI grids) as retrieved over Brussels by (a) TROPOMI and (b) MAX-DOAS. (c) Annual bias between tropospheric NO₂ VCD as observed by TROPOMI and MAX-DOAS (the negative values are shown with blue, zero with white, and positive with red). The black square shows the MAX-DOAS instrument location, the black polygon the national airport, the black dots the NO₂ hotspots, and the black line the Brussels Ring road.

Figure 22c shows the annual relative biases (e.g., $\mathrm{100}\times (\text{TROPOMI}-\text{MAX-DOAS})/\text{MAX-DOAS}$) per pixel. It is found that negative biases (i.e., MAX-DOAS larger than TROPOMI), ranging from −3 % to −38 %, are found for all the pixels.

The seasonal correlation plots for April 2018-February 2020 are displayed in Fig. 23. The highest correlation is found during spring (R=0.70), while lower correlations are reported in autumn, summer, and winter, with correlation coefficient values of 0.67, 0.61, and 0.57, respectively. It should be noted that during spring (2018 and 2019), the number of comparison points is smaller than for the other seasons because TROPOMI data start from the end of April 2018. During spring, the slope value is equal to 0.88, while during winter, summer, and autumn, the slope values are smaller (0.67, 0.45, and 0.56, respectively). Similar findings have been reported in several studies (Verhoelst et al., 2021; Tack et al., 2021; Judd et al., 2020; Dimitropoulou et al., 2020; Ialongo et al., 2019). When seasonally averaged TROPOMI and MAX-DOAS pixels (the pixels shown in Fig. 21) are compared one by one (see seasonal, SEAS, in Fig. 23), both the correlation coefficient (R in the range of 0.62–0.85) and slope values (s in the range of 0.36–1.84) improve considerably for spring and autumn.

Figure 23 Seasonal scatter plots of tropospheric NO₂ columns derived from the dual-scan MAX-DOAS and TROPOMI measurements over Brussels. Blue line: regression analysis results when all the MAX-DOAS and TROPOMI pixels are included in the comparison. Red line: seasonal average analysis generated by the pixels in Fig. 21.

The seasonal regression analysis parameters between TROPOMI and dual-scan MAX-DOAS measurements derived in the present study are compared to the same parameters presented in Dimitropoulou et al. (2020); see Table 4. Both studies make use of the dual-scan MAX-DOAS instrument in Uccle. In addition to the different approach (i.e., the retrieval of NO₂ horizontal profiles), in the present study, almost 2 years of measurements are used, while in Dimitropoulou et al. (2020), only 1 year is exploited for the TROPOMI validation. In Table 4, for the present study, only 1 year of measurements is used to have comparable time coverage for both studies. As presented in Table 4, the largest slope value is found in winter in both the present study and in Dimitropoulou et al. (2020). The season in which the highest correlation coefficient is obtained differs between the two studies (here in spring and in autumn in Dimitropoulou et al., 2020). The main advantage of the new mapping MAX-DOAS technique is the larger number of comparison points between TROPOMI and MAX-DOAS, leading to significantly more reliable statistics. In the present study, the deviation of the comparison points from the fitted regression line is increased, mainly because of the uncertainties in the horizontal inversion approach. The scatter increase is reflected in the correlation coefficient values, which are smaller for all seasons, except spring and winter. Regarding the slope value, it is larger during all seasons.

Table 4 Summary of the regression analysis parameters (e.g., correlation coefficient, R, and slope, s) and the number of data points (N) derived in the present study during only 1 year of observations (i.e., number of pixels) and in Dimitropoulou et al. (2020). Please note that R (seasonal) and s (seasonal) correspond to SEAS in Fig. 23.

7 Conclusions

A total of 2 years (March 2018 to February 2020) of dual-scan MAX-DOAS measurements in Uccle (urban background site located in the south of the Brussels-Capital Region) were used to develop a new strategy for the retrieval of near-surface NO₂ concentrations and aerosol extinction horizontal profiles. A full dual-scan measurement is composed of one vertical scan at a fixed azimuthal direction pointing towards the city center and horizontal scans in 10 azimuthal directions at a fixed low elevation angle (2^∘).

The first step of this new retrieval strategy is to analyze measured radiance spectra in six different fitting windows. This provides O₄ and NO₂ dSCDs at the following six wavelengths: 343, 360, 380, 447, 477, and 530 nm. Then, information about the vertical extent of NO₂ in the troposphere (MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$) is derived from profile retrievals in the main azimuthal direction performed using the OEM-based MMF algorithm. In the third step, a new parameterization technique is applied, with MLH${}_{{\mathrm{NO}}_{\mathrm{2}}}$, measured O₄ dSCDs, and measurement geometry being used as input parameters to retrieve the horizontal sensitivity of NO₂ and, consequently, the NO₂ near-surface concentrations, VCDs, near-surface aerosol extinction, and AODs in all the azimuthal directions for the six different wavelengths. Compared to the method presented in Dimitropoulou et al. (2020), the new retrieval method offers the possibility of the direct determination of $L_{{\mathrm{NO}}_{\mathrm{2}}}$ and near-surface aerosol extinction based on the measured O₄ dSCDs.

The retrieved dual-scan NO₂ near-surface concentrations and VCDs are verified via comparisons to the MMF NO₂ vertical profiles derived in the main azimuthal directions and in three additional azimuthal directions. Good overall agreement is found for the two comparisons during the 2 years of measurements.

The dependence of the horizontal sensitivity on the wavelength is then used to develop a new OEM-based horizontal distribution inversion approach. Considering a horizontal box model, horizontal NO₂ and aerosol extinction profiles are retrieved in an output horizontal grid of 500 m thickness starting from the instrument to each of the measurement maximum horizontal representative distances.

The daily variability of NO₂ horizontal profiles in all the azimuthal directions provides information about the location of the NO₂ hotspots in the Brussels-Capital Region and how the plumes are transported. Similarly, the NO₂ horizontal profiles' seasonal variability over March 2018–February 2020 reveals that the NO₂ hotspots are mainly found above the Brussels city center, the Drogenbos power plant, and the NW part of the Ring road during all seasons.

On 28 June 2019, airborne measurements (APEX) of NO₂ were performed over Brussels. The MAX-DOAS NO₂ VCD horizontal profiles are compared to APEX, car-mobile DOAS (i.e., AEROMOBIL), and TROPOMI measurements, and good overall agreement is found between the different datasets for this day.

In the last part of the study, MAX-DOAS retrievals are compared to TROPOMI tropospheric NO₂ observations over the March 2018–February 2020 period. The comparison of seasonal maps shows good overall agreement between the two datasets as to the NO₂ horizontal distribution over the Brussels area. Results also show that during all seasons, TROPOMI underestimates the MAX-DOAS tropospheric NO₂ columns.

Overall, our investigation about the spatial sampling led to the following two important findings.

1. The dual-scan multi-wavelength approach allows a good identification of the main emissions sources in urban regions, in agreement with the spatial allocation of the main emission sources observed by APEX and TROPOMI.

2. The characterization of the NO₂ concentration horizontal field using the dual scan multi-wavelength approach results in obtaining slope values closer to unity between TROPOMI and MAX-DOAS observations. The high spatial resolution of TROPOMI requires ground-based measurement that can provide information about the horizontal distribution of tropospheric NO₂ columns in urban regions.

To conclude, our study presents a new horizontal distribution inversion approach for NO₂ and aerosols developed by using dual-scan multi-wavelength MAX-DOAS measurements over an urban area. This approach provides a better characterization of the horizontal distribution of an important urban pollutant, NO₂, which leads to improved agreement between satellite and MAX-DOAS measurements in moderate to highly polluted conditions. Based on our study, further modifications of the measurement mode aiming at a better sampling of the vertical and horizontal NO₂ distribution could be implemented and investigated. For instance, performing vertical scans in several azimuthal directions throughout the day and/or horizontal scans in more than 10 azimuthal directions could further improve our knowledge about the tropospheric NO₂ spatial variability in urban regions and therefore the satellite validation results in those conditions.