Retrieval of tropospheric aerosol, NO2, and HCHO vertical profiles from MAX-DOAS observations over Thessaloniki, Greece: intercomparison and validation of two inversion algorithms

Abstract

In this study we focus on the retrieval of aerosol and trace gas vertical profiles from multi-axis differential optical absorption spectroscopy (MAX-DOAS) observations for the first time over Thessaloniki, Greece. We use two independent inversion algorithms for the profile retrievals: the Mexican MAX-DOAS Fit (MMF) and the Mainz Profile Algorithm (MAPA). The former is based on the optimal estimation method (OEM), while the latter follows a parameterization approach. We evaluate the performance of MMF and MAPA, and we validate their retrieved products with ancillary data measured by other co-located reference instruments. The trace gas differential slant column densities (dSCDs), simulated by the forward models, are in good agreement, except for HCHO, where larger scatter is observed due to the increased spectral noise of the measurements in the UV. We find an excellent agreement between the tropospheric column densities of NO${}_{\mathrm{2}}$ retrieved by MMF and MAPA ($\text{slope}=\mathrm{1.009}$, Pearson's correlation coefficient $R=\mathrm{0.982}$) and a good correlation for the case of HCHO ($R=\mathrm{0.927}$). For aerosols, we find better agreement for the aerosol optical depths (AODs) in the visible (i.e., at 477 nm) compared to the UV (at 360 nm), and we show that the agreement strongly depends on the O${}_{\mathrm{4}}$ scaling factor that is used in the analysis. The agreement for NO${}_{\mathrm{2}}$ and HCHO near-surface concentrations is similar to the comparison of the integrated columns with slightly decreased correlation coefficients. The seasonal mean vertical profiles that are retrieved by MMF and MAPA are intercompared, and the seasonal variation in all species along with possible sources is discussed. The AODs retrieved by the MAX-DOAS are validated by comparing them with AOD values measured by a CIMEL sun photometer and a Brewer spectrophotometer. Four different flagging schemes were applied to the data in order to evaluate their performance. Qualitatively, a generally good agreement is observed for both wavelengths, but we find a systematic bias from the CIMEL sun photometer and Brewer spectrophotometer measurements, due to the limited sensitivity of the MAX-DOAS in retrieving information at higher altitudes, especially in the UV. An in-depth validation of the aerosol vertical profiles retrieved by the MAX-DOAS is not possible since only in very few cases is the true aerosol profile known during the period of study. However, we examine four cases, where the MAX-DOAS provided a generally good estimation of the shape of the profiles retrieved by a co-located multi-wavelength lidar system. The NO${}_{\mathrm{2}}$ near-surface concentrations are validated against in situ observations, and the comparison of both MMF and MAPA revealed good agreement with correlation coefficients of $R=\mathrm{0.78}$ and $R=\mathrm{0.73}$, respectively. Finally, the effect of the O${}_{\mathrm{4}}$ scaling factor is investigated by intercomparing the integrated columns retrieved by the two algorithms and also by comparing the AODs derived by MAPA for different values of the scaling factor with AODs measured by the CIMEL sun photometer and the Brewer spectrophotometer.

1 Introduction

The planetary boundary layer (PBL), also called atmospheric boundary layer, is defined as the lowermost layer of the troposphere that is directly influenced by the terrestrial surface. The PBL height, at midlatitudes, expands typically up to 1–2 km during daytime , and its composition has a strong impact on weather, climate, and air quality. The increasing interest of understanding the PBL's structure and dynamics is apparent in various research fields, from air pollution analysis to weather prediction, and thus, continuous ground-based monitoring of both chemical composition and aerosol content of the PBL with high temporal resolution is of great importance.

Thessaloniki is a Mediterranean city, and it is the second largest city of Greece, located in the northern part of the country. Thessaloniki hosts approximately 10 % of the country's total population with more than 1 million inhabitants . With approximately 20 % of the country's industrial activity, it is considered one of the largest urban agglomerations in the Balkans . The air pollution sources in Thessaloniki are mainly industrial activities in the western part of the city as well as road transport and domestic heating during the cold period of the year, while the air quality of the city is affected by local topographic and meteorological characteristics . Nitrogen oxides (${\mathrm{NO}}_x=\mathrm{NO}+{\mathrm{NO}}_{\mathrm{2}}$), formaldehyde (HCHO) and aerosols are considered major atmospheric pollutants contained in the PBL of the city.

Nitrogen dioxide (NO${}_{\mathrm{2}}$) and HCHO are two important trace gas species of the atmosphere that play a critical role in tropospheric photochemistry , participating in the formation of tropospheric ozone (O${}_{\mathrm{3}}$), while aerosols can have a strong influence on air quality and climate through effects on radiation . Both NO${}_{\mathrm{2}}$ and HCHO are toxic to humans in high concentrations and can lead to severe health conditions. HCHO is a short-lived product derived by the oxidation of volatile organic compounds (VOCs). Its sources are both natural (i.e., oxidation of VOCs emitted from plants) and anthropogenic (i.e., biomass burning, industry-related emissions, and road transport) . NO${}_{\mathrm{2}}$ is mainly produced by the oxidation of nitrogen monoxide (NO), and in most urban areas its sources include fossil fuel combustion, biomass burning, soil emissions, and lightning . Moreover, under certain meteorological conditions, NO${}_{\mathrm{2}}$ may participate in the formation of secondary aerosols . Given the influence of NO${}_{\mathrm{2}}$, HCHO, and aerosols on air quality and climate, it is of high environmental and research importance to accurately and continuously monitor their spatio-temporal distribution in the troposphere.

Multi-axis differential optical absorption spectroscopy (MAX-DOAS) is a well-established ground-based passive remote sensing technique that received considerable attention during the past decades and is nowadays widely used in many studies in order to simultaneously detect trace gases and aerosols mainly in the PBL and in the lowermost free troposphere (e.g., , and references therein). Such trace gases include NO${}_{\mathrm{2}}$, HCHO, sulfur dioxide (SO${}_{\mathrm{2}}$), water vapor (H${}_{\mathrm{2}}$O), ozone (O${}_{\mathrm{3}}$), nitrous acid (HONO), iodine oxide (IO), glyoxal (CHOCHO), and bromine oxide (BrO). The MAX-DOAS measurement technique utilizes scattered sunlight in the ultraviolet (UV) and visible (VIS) parts of the electromagnetic spectrum received from different elevation angles, and the measured spectra are analyzed by differential optical absorption spectroscopy (DOAS) for the determination of the differential slant column densities (dSCDs). Information about the vertical distribution of aerosols and trace gases can be retrieved from a single elevation sequence (i.e., spectra recorded at different elevation angles that belong to the same azimuthal direction) using suitable inversion algorithms. The products retrieved by the inversion algorithms include, among others, estimates of the profile shape, tropospheric vertical column densities (VCDs), and near-surface concentrations.

Nowadays, there is a variety of such inversion algorithms for the retrieval of vertical profiles from MAX-DOAS measurements using different techniques. These algorithms are mainly separated into those that retrieve the profiles based on the optimal estimation method (OEM) and into those that rely on a few parameters to characterize the atmospheric profile (parameterization approach). Both OEM-based and parameterized inversion algorithms have been tested and intercompared so far in many studies using either synthetic data e.g., or actual MAX-DOAS measurements, as for example, during the Cabauw Intercomparison of Nitrogen Dioxide Measuring Instruments 2 (CINDI-2) campaign . Here, we use two of the already tested inversion algorithms to analyze MAX-DOAS measurements conducted at Thessaloniki, Greece, for the retrieval of aerosol, NO${}_{\mathrm{2}}$, and HCHO vertical profiles and column densities. These algorithms are the Mexican MAX-DOAS Fit (MMF) v2020_04 and the Mainz Profile Algorithm (MAPA) v0.98 . The former is based on the OEM, while the latter follows the parameterization approach, and both are adopted by the Fiducial Reference Measurements for Ground-Based DOAS Air-Quality Observations (FRM${}_{\mathrm{4}}$DOAS) project (https://frm4doas.aeronomie.be/, last access: 5 March 2021). In this work we evaluate the performance of the two algorithms, we validate their results with reference datasets, and we investigate the effect of applying different flagging schemes to the retrieved products. Additionally, by using two independent inversion algorithms, we aim at producing a reference MAX-DOAS dataset of higher quality for further research activities in Thessaloniki (e.g., validation of satellite-retrieved tropospheric products). Thessaloniki is also part of the FRM${}_{\mathrm{4}}$DOAS project, which aims at the development of the first central processing system for MAX-DOAS observations. Even though the measured spectra are regularly submitted and analyzed on a near-real-time basis, in this work both MMF and MAPA runs are performed offline in order to obtain more flexibility in the analysis and also to investigate and optimize the retrieval settings, particularly for Thessaloniki.

The article is structured as follows. In Sect.  the instrumentation, the MAX-DOAS retrieval settings, and a brief description of the profiling algorithms are reported, along with the methodology used in this analysis. In Sect.  we present the results of the comparison between different products retrieved by MMF and MAPA. In Sect.  the validation results of the retrieved products with ancillary data are presented, and in Sect.  the main conclusions of this article are summarized.

2 Data and methodology

2.1 Instrumentation

A 2D MAX-DOAS system (Phaethon) operates regularly on the rooftop (20 m above ground) of the Physics Department building of the Aristotle University of Thessaloniki (40.634${}^{\circ }$ N, 22.956${}^{\circ }$ E), about 60 m above sea level (a.s.l.). The measurement site is located near the city center of Thessaloniki (Fig. ). The prototype system was developed in 2006 at the Laboratory of Atmospheric Physics (LAP) and has been upgraded ever since for the retrieval of tropospheric NO${}_{\mathrm{2}}$ VCDs and total ozone columns . The current version of the system comprises a single-channel ultra-low stray-light AvaSpec-ULS2048x64-EVO ($f=\mathrm{75}$ mm) spectrometer by Avantes, the entrance optics, and a two-axis tracker. The spectrometer's detector is a back-thinned Hamamatsu charge-coupled device (CCD) array of 2048 pixels with signal-to-noise ratio (SNR) $\mathrm{450}:\mathrm{1}$ for a single measurement at full signal. The spectrometer covers the spectral range 280–539 nm and uses a 50 $\mathrm{\mu }\mathrm{m}$ wide entrance slit. Mercury discharge lamp spectra were recorded to determine the instrument's slit function, and the spectral resolution was found to be $\sim \mathrm{0.55}$ nm full width at half maximum (FWHM) at 436 nm. The spectrometer is positioned inside a thermally isolated box, where the temperature is maintained at $+\mathrm{10}$ ${}^{\circ }$C using a thermoelectric Peltier system. The entrance optics are mounted on a two-axis tracker with two stepper motors controlling the azimuth viewing angle ($\mathrm{0}{}^{\circ }\le \mathit{\varphi }\le \mathrm{360}{}^{\circ }$) and the elevation viewing angle ($\mathrm{0}{}^{\circ }\le \mathit{\alpha }\le \mathrm{90}{}^{\circ }$) with pointing resolution of 0.125${}^{\circ }$, allowing both direct-sun and off-axis observations. A third motor rotates a filter wheel of eight positions with different optical components (diffuser, attenuation, and band-pass filters), used for the measurement of direct-sun and scattered radiation spectra and an opaque position for the measurement of the dark signal. The instrument operates automatically and is controlled by a custom-made software, developed at LAP. The entrance also comprises a telescope with a plano-convex lens that focuses the collected solar radiation onto one end of an optical fiber. The system's field of view (FOV) was characterized using a distant light source and was found $\sim \mathrm{1}$${}^{\circ }$. Simultaneous azimuth and elevation angle calibration are regularly performed by sighting the sun, so no horizon scans are necessary for the elevation angle calibration.

Figure 1

The Phaethon MAX-DOAS system in the middle and a panoramic view (east–south–west) of the measurement site.

[Figure omitted. See PDF]

A routine MAX-DOAS measurement cycle starts by orienting the optics at a certain azimuth viewing direction followed by the measurement of scattered radiation spectra at the elevation angles: 90 (zenith), 30, 15, 12, 10, 8, 6, 5, 4, 3, 2, and 1${}^{\circ }$ in this order. For this study, the system was configured to measure at four consecutive azimuth angles of 142, 185, 220, and 255${}^{\circ }$, illustrated in Fig.  with arrows of different colors. Based on the intensity of the measured spectra during an elevation scan at 142${}^{\circ }$ azimuth, the viewing direction of the 1${}^{\circ }$ elevation angle was found to be partly blocked by obstacles, such as trees and buildings in the campus. Thus, $\mathit{\alpha }=\mathrm{1}$${}^{\circ }$ in this particular direction was excluded from the profiling analysis. In order to achieve high SNR values and to avoid saturated spectra, the number of scans of each individual measurement and the exposure time of the CCD are automatically adjusted by the operating software according to the received intensity by the detector. The integration time at each elevation angle is $\sim \mathrm{60}$ s, and a full measurement sequence for all azimuth directions lasts about 1 h.

Figure 2

Location of the MAX-DOAS system (white dot) and the in situ NO${}_{\mathrm{2}}$ measurement site (yellow dot). The arrows in different colors represent the azimuth viewing directions, $\mathit{\varphi }$, of the MAX-DOAS observations (i.e., purple: 142${}^{\circ }$; blue: 185${}^{\circ }$; green: 220${}^{\circ }$; red: 255${}^{\circ }$). The base map is taken from © Google Maps, https://www.google.com/maps/@40.6194132,22.9538922,16768m/data=!3m1!1e3!5m1!1e4 (last access: 3 March 2022).

[Figure omitted. See PDF]

2.2 MAX-DOAS measurements and slant column retrieval settings

The primary retrieved product from the analysis of the measured MAX-DOAS spectra is the dSCD of several trace gases at different elevation angles. The dSCD of a trace gas at an elevation angle $\mathit{\alpha }$ (${\text{dSCD}}_{\mathit{\alpha }}$) can be calculated as the difference between the slant column density, i.e., its concentration integrated along the light path (${\text{SCD}}_{\mathit{\alpha }}$) and the SCD of a Fraunhofer reference spectrum (FRS), usually measured at the zenith (${\text{SCD}}_{\text{ref}}$):

1$${\text{dSCD}}_{\mathit{\alpha }}={\text{SCD}}_{\mathit{\alpha }}-{\text{SCD}}_{\text{ref}}.$$

The MAX-DOAS spectra that are used in this study have been recorded for 1 year (from May 2020 through May 2021), and a zenith spectrum is selected as the FRS in order to account for the Fraunhofer lines and the stratospheric contribution of the absorbers . Since the system is scheduled to perform both direct-sun and MAX-DOAS observations during the day, the zenith spectra of two consecutive elevation sequences may have a large time difference (duration of the first sequence plus the duration of two direct-sun measurements). So, in this study, the zenith spectrum of each sequence was selected as the FRS for the DOAS-based retrieval of the collision-induced oxygen complex (O${}_{\mathrm{2}}$–O${}_{\mathrm{2}}$ or O${}_{\mathrm{4}}$) and the trace gas dSCDs and not the average or the time-interpolated spectrum between the zenith spectra of the two consecutive sequences. The dSCDs of O${}_{\mathrm{4}}$ and trace gases are derived from the recorded spectra by applying the DOAS technique , while the measured spectra are analyzed using the QDOAS (version 3.2, September 2017) spectral fitting software suite developed by BIRA-IASB (https://uv-vis.aeronomie.be/software/QDOAS/, last access: 5 March 2021) . The retrieval settings are based on results from the CINDI-2 campaign (http://www.tropomi.eu/data-products/cindi-2/, last access: 5 March 2021) e.g.,, the Quality Assurance for Essential Climate Variables (QA4ECV) project (http://www.qa4ecv.eu/, last access: 5 March 2021), and the Network for Detection of Atmospheric Composition Change (NDACC) protocol for UV–VIS measurements (http://www.ndaccdemo.org/data/protocols/, last access: 5 March 2021). The spectral retrieval settings and the trace gas absorption cross sections that are included in the DOAS fit are listed in Table . The wavelength calibration of the measured spectra is achieved by shifting and stretching them against a highly resolved solar reference spectrum . Even though the spectrometer is operating in a temperature-controlled environment, small diurnal temperature variations may occur. Thus, dark spectra are measured after each elevation sequence for all of the exposure times that were used during the sequence. This procedure might be time-consuming, but it assures that the solar and dark spectra are measured at the same temperature. The dark spectra are then subtracted from the scattered radiation spectra prior to the DOAS analysis. Figure  shows a typical example of the DOAS analysis of a spectrum recorded on 9 July 2020 at 07:50 UTC at 3${}^{\circ }$ elevation angle. During the whole period of study no apparent system-related issue or instrument degradation is observed.

Figure 3

A typical example of the DOAS retrieval of NO${}_{\mathrm{2}}$, HCHO, O${}_{\mathrm{4}}$ (VIS), and O${}_{\mathrm{4}}$ (UV) dSCDs derived from a MAX-DOAS measurement on 9 July 2020 at 07:50 UTC ($\text{SZA}=\mathrm{38.85}{}^{\circ }$) at 3${}^{\circ }$ viewing elevation angle. The DOAS fits are presented in the subpanels of panel (a). The black lines represent the measured spectra, and red lines are the fitted O${}_{\mathrm{4}}$ and trace gas cross sections. The subpanels of panel (b) show the residual of the DOAS fits.

[Figure omitted. See PDF]

2.3 Retrieval of the vertical profile

The retrieval of vertical profiles (extinction and concentration profiles for aerosols and trace gases, respectively) from MAX-DOAS measurements typically involves three major steps , independent of the retrieval approach. In the first step, the O${}_{\mathrm{4}}$ dSCDs and the trace gas dSCDs (in this case NO${}_{\mathrm{2}}$ and HCHO) are derived by applying the DOAS fitting technique to the measured spectra, as described in Sect. . Next, the O${}_{\mathrm{4}}$ dSCDs retrieved for each elevation angle of the same sequence are used as input to the algorithm for the retrieval of the aerosol extinction vertical profile. In the end, the trace gas dSCDs are used as input to the algorithm for the retrieval of the trace gas vertical profile, along with the aerosol extinction profile, calculated in the previous step.

As mentioned already, the profiling algorithms that have been developed so far and are commonly used within the MAX-DOAS community are either based on the OEM or follow the parameterization approach. How the O${}_{\mathrm{4}}$ and trace gas dSCDs are handled for the retrieval of the vertical profiles depends on each algorithm's approach. However, the principal idea of both OEM and parameterized inversion algorithms is the same: a layered model atmosphere with defined parameters is assumed in a forward radiative transfer model (RTM), and it is used in order to simulate the O${}_{\mathrm{4}}$ and trace gas dSCDs, taking into account the viewing geometry, i.e., the solar zenith angle (SZA), the elevation angle, and the relative azimuth angle. The forward models and how the dSCDs are simulated are described in for MAPA and in for MMF. The extinction and concentration vertical profiles are derived by inverting the forward model, i.e., by finding the model parameters, for which the difference between the simulated and the measured dSCDs is minimized, based on a cost function.

Table 1

DOAS fit settings for NO${}_{\mathrm{2}}$, HCHO, O${}_{\mathrm{4}}$ (VIS), and O${}_{\mathrm{4}}$ (UV).

Parameter	Data source	Trace gas
		NO${}_{\mathrm{2}}$ and O${}_{\mathrm{4}}$	HCHO	O${}_{\mathrm{4}}$
		(VIS)		(UV)
Spectral range		425–490 nm	324.5–359 nm	338–370 nm
NO${}_{\mathrm{2}}$ (298 K)	, $I_{\mathrm{0}}$-corrected	✓	✓	✓
	(SCD in ${\mathrm{10}}^{\mathrm{17}}$ molec. cm${}^{-\mathrm{2}}$)
NO${}_{\mathrm{2}}$ (220 K)	, $I_{\mathrm{0}}$-corrected	✓		✓
	(SCD in ${\mathrm{10}}^{\mathrm{17}}$ molec. cm${}^{-\mathrm{2}}$)
O${}_{\mathrm{3}}$ (223 K)	, $I_{\mathrm{0}}$-corrected	✓	✓	✓
	(SCD in ${\mathrm{10}}^{\mathrm{20}}$ molec. cm${}^{-\mathrm{2}}$)
O${}_{\mathrm{3}}$ (243 K)	, $I_{\mathrm{0}}$-corrected		✓	✓
	(SCD in ${\mathrm{10}}^{\mathrm{20}}$ molec. cm${}^{-\mathrm{2}}$)
O${}_{\mathrm{4}}$ (293 K)		✓	✓	✓
BrO (223 K)			✓	✓
HCHO (297 K)			✓	✓
H${}_{\mathrm{2}}$O (296 K)	HITEMP,	✓
Ring	Ring spectra calculated by QDOAS according to	✓	✓	✓
Polynomial degree		5	5	5
Intensity offset		Constant	Order 1	Constant
Wavelength calibration	Based on a high-resolution solar reference spectrum

2.4 MMF

The Mexican MAX-DOAS Fit (MMF) v2020_04 is an OEM-based profiling algorithm that relies on online RTM simulations using VLIDORT version 2.7 as forward model. The input parameters for each atmospheric layer are calculated from temperature and pressure profiles, the trace gas concentration in each layer, and the aerosol properties. The aerosol properties, which are the same for all layers, are the single-scattering albedo (SSA) and the asymmetry parameter (using the Henyey–Greenstein phase function, , to calculate the phase function moments). Furthermore, the wavelength of the retrieval and the surface albedo need to be specified as additional input parameters. The retrieval algorithm comprises an aerosol extinction profile retrieval and a trace gas profile retrieval. The former constrains the aerosol extinction profile in the forward model of the trace gas retrieval. The inversion uses constrained damped least-square fitting with an optimal estimation regularization. In the used version, both the a priori and the covariance matrix are constructed. More details about the a priori settings and the input parameters can be found in Sect. . The retrieval algorithm provides the aerosol extinction profiles, trace gas partial column profiles, their integrated quantities, the corresponding noise and smoothing errors, the averaging kernel, the degrees of freedom, and a quality flag of the retrieval. The quality flagging system of MMF is based on the convergence of the algorithm, the root mean square of the difference between measured and simulated dSCDs, the reported degrees of freedom, and the stability of the retrieval.

2.5 MAPA

The Mainz Profile Algorithm (MAPA) v0.98 is a profiling algorithm developed by the Max Planck Institute for Chemistry (MPIC) that is based on a parameterization approach. MAPA does not rely on online RTM simulations, but its forward model is provided as pre-calculated differential air mass factor (dAMF) look-up tables (LUTs) at multiple wavelengths. These LUTs have been calculated offline by a full spherical RTM, McArtim , following a backward Monte Carlo approach. Just like MMF, MAPA is based on a two-step process in order to retrieve the aerosol and trace gas vertical profiles. It uses three main parameters to characterize the atmospheric profile: the column parameter, $c$ (i.e., AOD for aerosols and VCD for trace gases); the layer height, $h$; and the shape parameter, $s$. Additionally, a fourth optional parameter can be included, the O${}_{\mathrm{4}}$ scaling factor, which was initially introduced by in order to achieve agreement between the measured dSCDs and the forward model simulations. Unlike MMF, MAPA is not based on the OEM, so no a priori assumption of the vertical profile is required. In some cases this can be an advantage since a priori information and constraints are usually difficult to estimate. MAPA also provides a detailed flagging algorithm that is based on thresholding techniques applied to different parameters, in order to evaluate whether the retrieved profile can be trusted. By default (and within this study), the flags are identical for the species retrieved in the UV and VIS spectral range. The flags that are defined in MAPA v0.98 are mainly based on the agreement between the measured and modeled dSCDs, the consistency of the derived Monte Carlo parameters, and the shape of the profile. More details about MAPA and its flagging algorithm can be found in .

2.6 Input parameters and settings

During MAPA calculations, depending on the aerosol or trace gas retrieval, a LUT corresponding to the central wavelength of the O${}_{\mathrm{4}}$ or trace gas fitting window is selected (i.e., 360 nm for O${}_{\mathrm{4}}$ in the UV, 343 nm for HCHO, 460 nm for NO${}_{\mathrm{2}}$, and 477 nm for O${}_{\mathrm{4}}$ in the VIS). These wavelengths are also used in the RTM simulations of MMF. For the calculation of the dAMF LUTs, MAPA's radiative transfer simulations were performed with a typical fixed set of parameters for all wavelengths , which can describe the majority of all potential measurement sites. MMF, on the other hand, relies on online RTM simulations, and so the aerosol and surface parameters can be adjusted to the most suitable values. In this study, the aerosol optical properties that are used as input for the simulations of MMF are based on 15 years of climatological data measured by a co-located CIMEL sun photometer. Figure  shows the frequency distribution of the Ångström exponent, AOD, asymmetry factor, and SSA in Thessaloniki, while their values that are used as input to each inversion algorithm are listed in Table . Discrepancies between MMF and MAPA due to small differences in these selected parameters are expected to be minor.

Table 2

The RTM settings that were used in MMF and MAPA for Thessaloniki.

Parameter	Inversion algorithm
	MAPA	MMF
Aerosol single-scattering albedo	0.95	0.95
Aerosol asymmetry parameter	0.68	0.69
Surface albedo	0.05	0.06
Ångström exponent	1.4	1.4

MMF requires a priori profile and covariance matrix information for the profile retrievals. The “a priori” term represents knowledge of the true state before the measurement is performed. However, the true shape of the trace gas vertical profiles at Thessaloniki is generally not known, while the true state of the aerosol profiles is known only in certain cases during the period of study. Thus, the retrieval is based on constructed exponentially decreasing a priori profiles with a scale height of 1 km, which are considered a reasonable estimate of the true profiles. Since no covariance matrix information is available, the covariance matrix is also constructed from the a priori profile. The AOD as well as the trace gas VCDs in Thessaloniki vary substantially throughout the year. In order to take into account the annual variability, we use the square of 50 % of the a priori on the diagonal elements of the covariance matrix for aerosols and 100 % for NO${}_{\mathrm{2}}$ and HCHO. The loose constraint of the latter is due to the higher variability of the trace gas vertical columns over the course of the year. Both for aerosols and trace gases, the off-axis elements of the covariance matrix were constructed by assuming a Gaussian function with a correlation length of 200 m, as described in . Additionally, based on empirical tests, the progress of the convergence is faster when using an a priori VCD or AOD below the true value for reasons that are not yet identified. Thus, the a priori AODs were set to 0.25 and 0.15 for the aerosol retrievals at 360 and 477 nm, respectively. For the trace gas retrievals we have used a priori VCDs of $\mathrm{4}\times {\mathrm{10}}^{\mathrm{15}}$ and $\mathrm{6}\times {\mathrm{10}}^{\mathrm{15}}$ molec. cm${}^{-\mathrm{2}}$ for NO${}_{\mathrm{2}}$ and HCHO, respectively, based on data derived from the MAX-DOAS by applying the geometrical approximation to the dSCDs measured at 30 and 15${}^{\circ }$ elevation angles. The LUTs used in MAPA cover the following ranges: 0–5 for the AOD, 0.02–5 km for the layer height, and 0.2–1.8 for the profile shape parameter .

The temperature and pressure vertical profiles that are used as input in this study are identical for both MMF and MAPA. We have used climatological profiles for Thessaloniki produced by MPIC that are based on $\sim \mathrm{16}$ years of re-analysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF). The temperature and pressure profiles are interpolated to the day and time of each elevation sequence. We have also tried to use temperature and pressure profiles measured by radiosondes, launched on a daily basis at Thessaloniki Airport ($\sim \mathrm{13}$ km away from the measurement site), as input, but since no major effect is observed on the retrieved products, these results are not presented. Both algorithms are configured to export the retrieved vertical profiles to the same output grid ranging from the ground up to 4 km with 200 m vertical resolution.

As already mentioned, the recorded spectra are also analyzed by a central processing system in the frame of the FRM${}_{\mathrm{4}}$DOAS project. The analysis is carried out using default values of several parameters, which are reasonable for all potential measurement sites, while in this study we try to optimize the performance of MMF and MAPA in particular for Thessaloniki (see discussion above). In the FRM${}_{\mathrm{4}}$DOAS analysis a time-interpolated spectrum between the zenith spectra of two consecutive elevation sequences is used as the FRS for the dSCDs retrievals. Thus, the dSCDs that are used for the retrieval of the vertical profiles are slightly different than those used in the current study. The default FRM${}_{\mathrm{4}}$DOAS settings include SSA of 0.92 and an asymmetry factor of 0.68. Yet, such small differences should have a negligible effect on the retrieved vertical profiles. The Ångström exponent is set to 1, and the same a priori aerosol extinction vertical profile (AOD of 0.18) is used for the retrievals in both the UV and VIS spectral ranges. The covariance matrices are constructed from the a priori profiles, but the square of 50 % of the a priori is used on the diagonal elements of the covariance matrices for all species. Currently, the partly blocked elevation angle of 1${}^{\circ }$ at 142${}^{\circ }$ azimuth (Sect. ) is not excluded from the analysis. MAPA retrievals are performed using three different O${}_{\mathrm{4}}$ scaling factors (i.e., 0.8, 1.0, and a variable scaling factor). In order to further investigate the effect of the O${}_{\mathrm{4}}$ scaling factor (see Appendix ) in this study, we include an extra value of 0.9.

Figure 4

Frequency distribution of the Ångström exponent, AOD, asymmetry factor, and SSA measured by a CIMEL sun photometer in Thessaloniki for the period 2005–2021.

[Figure omitted. See PDF]

2.7 Ancillary data

This section briefly describes the supporting instruments that are used in this study for comparison and validation of MAX-DOAS-derived products. The ancillary data include measurements of a CIMEL sun photometer, a Brewer spectrophotometer, an aerosol lidar system, and an in situ NO${}_{\mathrm{2}}$ monitoring station. Except for the in situ NO${}_{\mathrm{2}}$, all remote sensing instruments that are used in this study (i.e., the MAX-DOAS, CIMEL sun photometer, Brewer spectrophotometer, and lidar) are located at the same measurement site, about 60 m a.s.l. The effect of the different viewing geometries and the retrieval techniques that each system utilizes are discussed in the corresponding following sections.

2.7.1 CIMEL sun photometer

Since 2003 a sun–sky photometer (CIMEL) has provided spectral measurements of the AOD at Thessaloniki as part of the NASA's Aerosol Robotic Network (AERONET) (https://aeronet.gsfc.nasa.gov/, last access: 5 March 2021). The CIMEL sun photometer is an automated, well-calibrated scanning filter radiometer specifically developed for the retrieval of the AOD at seven wavelengths (i.e., 340, 380, 440, 500, 670, 870, and 1020 nm) by using direct-sun observations. The technical specifications of the instrument are given in . The instrument is calibrated regularly following the procedures and the guidelines of AERONET. The AERONET database provides three distinct levels for data quality. Level 1.0 is defined as pre-screened data (i.e., no quality assurance criteria are applied). Version 3 of Level 1.5 represents near-real-time automatic cloud screened data, while Level 2.0 applies additional pre- and post-field calibrations. In this paper, we use the AERONET Level 1.5 data, since the Level 2.0 data for the period of study are not yet published. In order to compare with the AOD retrieved by the MAX-DOAS, the AODs at 360 and 477 nm have been calculated using the Ångström exponent between the standard spectral bands of the instrument.

2.7.2 Brewer spectrophotometer

The Brewer spectrophotometer with serial number 086 (B086) is a double monochromator that has performed spectrally resolved measurements of the direct and global solar irradiance at Thessaloniki since 1993 . The wavelength range of B086 is 290–365 nm, and its spectral resolution is 0.55 nm at full width at half maximum (FWHM). The wavelength calibration is performed by scanning the emission lines of spectral discharge lamps, while maintenance of the absolute calibration is achieved by regularly scanning the spectral irradiance of a calibrated 1000 W quartz–halogen tungsten lamp .

Although the Brewer spectrophotometer's initial purpose was the retrieval of total ozone columns, research activities have shown that the spectral AOD can be calculated from direct irradiance measurements by following two main approaches: the first is based on the absolute calibration of the direct-sun spectra measured by the Brewer spectrophotometer , while the second uses the Langley extrapolation method (relative calibration) . In both cases the spectral AOD is calculated as the residual optical depth after subtracting the optical depths due to molecular scattering and the O${}_{\mathrm{3}}$ and SO${}_{\mathrm{2}}$ absorption from the total atmospheric optical depth . Since 1997, the direct solar irradiance spectra measured by the Brewer spectrophotometer have been calibrated , so in this study we use the former approach (i.e., absolute calibration) for the retrieval of the spectral AOD. In order to compare with the AOD retrieved by the MAX-DOAS, the AOD at 477 nm is calculated using climatological monthly mean values of the extinction Ångström exponent derived from measurements of the CIMEL sun photometer in Thessaloniki. Details on the procedure of the Brewer spectrophotometer's direct solar irradiance spectra absolute calibration as well as the spectral AOD retrieval methodology can be found in , , and .

2.7.3 Lidar

Thessaloniki has been a member station of the European Lidar Aerosol Network (EARLINET, https://www.earlinet.org, last access: 5 March 2021) since 2000, providing regular aerosol profile measurements, following EARLINET’s schedule (Monday morning and Monday and Thursday evening), during extreme events and at satellite overpasses (e.g., AEOLUS, OMI).

THEssaloniki LIdar SYStem (THELISYS) is a multi-wavelength Raman/depolarization lidar system, which has been gradually upgraded regarding its operational wavelengths and the detection configuration. All the quality standards, established within EARLINET, are followed in order to assure the high quality of the THELISYS products, which are publicly available in the EARLINET database (https://www.earlinet.org/index.php?id=125, last access: 5 March 2021). A detailed description of THELISYS technical specifications and algorithm can be found in .

The final products derived from the raw lidar data processing are the aerosol backscatter coefficient at 355, 532, and 1064 nm; the aerosol extinction coefficient at 355 and 532 nm; and the linear particle $/$ volume depolarization ratio at 532 nm. During the day, the data acquisition is limited to the signals that arise from the elastic scattering of the laser beam by the air molecules and the atmospheric aerosol. The Klett–Fernald algorithm in backward integration mode is applied , and the backscatter coefficient profiles are produced. Constant a priori climatological values of the ratio between the extinction and the backscatter coefficient (lidar ratio) were assumed in this daytime method. Values of 60, 50, and 40 were used for 355, 532, and 1064 nm, respectively, given the atmospheric situations that occur over Thessaloniki . The resulting uncertainties are discussed in depth by and can be as high as 50 % if there is no information about the actual lidar ratio, during extreme atmospheric conditions.

Another source of uncertainty during the lidar signal processing is the system's overlap function, which determines the altitude, above which a profile contains trustworthy values. In our analysis, the correction is not available for the daytime retrievals. Thus, an overlap function from the previous nighttime measurement or a mean overlap profile is applied. The starting height is set to the full overlap height (approximately 0.6 km), assuming height-independent backscatter below 0.6 km, equal to the backscatter measured at this height, to account for both the incomplete overlap within the lidar profile and atmospheric variability in the lowermost tropospheric part. This overlap effect generally introduces uncertainties in the calculation of the columnar products (e.g., AOD). However, long-term comparisons have shown similar decreasing trends of the AOD at 355 nm between the EARLINET and the AERONET datasets ($-\mathrm{23.2}$ % per decade and $-\mathrm{22.3}$ % per decade, respectively). The AODs at 355 nm measured by the lidar have also been compared with the Brewer spectrophotometer's retrievals, showing a generally good correlation of 0.7 .

2.7.4 In situ

Near-surface concentrations of different air pollutants, including NO${}_{\mathrm{2}}$, NO, SO${}_{\mathrm{2}}$, CO, and O${}_{\mathrm{3}}$, are measured in Thessaloniki by in situ instruments as part of the Network for Air Quality Monitoring of the Municipality of Thessaloniki. NO${}_{\mathrm{2}}$ is being monitored by chemiluminescence detectors that are mainly distributed around the city center. Most of the network stations are installed very close to the ground (sampling inlet at $\sim \mathrm{3}$ m) and are strongly affected by local traffic emissions. In this study, we use hourly mean (which is the highest available temporal resolution) in situ NO${}_{\mathrm{2}}$ concentrations measured at the “Eftapyrgion” site (40.644${}^{\circ }$ N, 22.957${}^{\circ }$ E, 174 m a.s.l.), which is located in an urban background area at a distance of $\sim \mathrm{1.2}$ km from the MAX-DOAS system to the north (Fig. ). The in situ measurements, spanning from May 2020 to March 2021, are used in order to validate the MAX-DOAS-derived NO${}_{\mathrm{2}}$ near-surface concentrations. Even though this site is located opposite to the MAX-DOAS system's azimuth viewing directions, it has been selected because the vertical and horizontal displacement of the two instruments is small, but also because it is the only site of the network almost unaffected by local traffic emissions and therefore can be considered more representative of the average NO${}_{\mathrm{2}}$ concentrations in the local boundary layer.

3 Results and discussion

In this section, we present results of the trace gas and aerosol quantities retrieved by the two inversion algorithms. We intercompare the dSCDs simulated by the forward models, the integrated columns (i.e., VCDs and AODs for trace gases and aerosols, respectively), the near-surface concentrations, and the seasonal mean vertical profiles between MMF and MAPA. Since MAPA is based on a parameterization approach, no information about averaging kernels is provided; hence, results on averaging kernels are presented only for MMF.

The MAX-DOAS system operates at a site where the northern viewing directions are blocked by buildings of the campus and the city, so the system is configured to perform sequences of elevation scans at azimuth directions in the southern sector, as illustrated in Fig. . As a result, scattered radiation spectra may be measured during the day at azimuths close to the solar azimuth angle. In such cases, RTM simulations might face difficulties in properly calculating the dAMF due to increased aerosol forward scattering, usually leading to underestimation of the true dAMF. For small scattering angles the uncertainties caused by the incorrect description of the phase function can also become important, and the results for such viewing geometries should be treated with caution. Therefore, the elevation sequences measured at azimuth angles relative to the sun of less than 5${}^{\circ }$ are excluded from the analysis. In addition, the elevation sequences, for which the retrieved AOD from the MAX-DOAS inversion algorithms is greater than 1.5, are filtered out, since such high aerosol loads are unrealistic for Thessaloniki (Fig. ). Negative columns can occur in the trace gas retrievals of MAPA within the Monte Carlo ensemble, and they are by default not removed, but this is not possible for MMF retrievals since, in its current version, MMF operates in logarithmic state vector space. For NO${}_{\mathrm{2}}$, no valid negative columns are retrieved, but for HCHO, MAPA reports negative columns for $\sim \mathrm{8.5}$ % of the valid data. In order to compare meaningful results between the two algorithms, the negative columns are removed from the initial dataset.

The individual flagging schemes of MMF and MAPA have been discussed elsewhere. Based on synthetic data, reported that the quality flagging criteria of MAPA might be too strict, since a large fraction of data were flagged as invalid, even though the algorithm successfully removed almost all outliers. In our study, MAPA flags a larger fraction of data as invalid, compared to MMF, for all the retrieved species. The percentage of the valid data flagged by MAPA and MMF (individually and combined) is presented in Table . Since in MAPA retrievals no a priori constraints are used, more strict flagging needs to be applied for retrieved dSCDs that are characterized by large uncertainties (e.g., due to larger fit error or the effects of clouds). Especially for HCHO, the apparent worse performance of MAPA could be explained by the lower SNR in the UV, along with the higher HCHO profile height compared to NO${}_{\mathrm{2}}$ (see discussion in Sect. ) and the decreasing sensitivity towards higher altitudes. The retrieval results are sensitive to the validity flagging approach, which is further investigated in the next section. No cloud filtering is applied to the data prior to the profiling analysis. Neither MMF nor MAPA include a direct cloud flagging system. However, some flags that are included in the flagging algorithms of MMF and MAPA are sensitive to clouds. Hence, in order to achieve retrievals of high quality and to ensure that the MAX-DOAS measurements performed under broken cloud conditions are filtered out, an elevation sequence is considered valid as long as it is flagged as valid by both MMF and MAPA. This is the default flagging scheme for NO${}_{\mathrm{2}}$, HCHO, and AOD at 477 nm, and all the results shown in the next sections follow this flagging approach unless stated otherwise. For AOD at 360 nm the flags reported by MAPA are considered default, since this approach performs better when comparing the MAX-DOAS results with other reference instruments (Sect. ), although the reason for this behavior has not yet been identified. Also, since the issue for selecting the optimum O${}_{\mathrm{4}}$ scaling factor remains unresolved , we let MAPA determine an optimum O${}_{\mathrm{4}}$ scaling factor (variable) for each elevation sequence, and this option is selected as the default for the retrievals.

It should be noted that in the following sections an orthogonal distance regression (ODR or bivariate least-squares) has been used instead of an ordinary linear regression (OLR or standard least-squares) for the comparison of the retrieved products derived by MMF and MAPA, in order to equally treat the two algorithms since none of them depend on the other. The discrepancies in the regression slopes and intercepts arising in the OLR when comparing independent variables, and the appropriateness of ODR, are discussed in . The ODR results are also sensitive to the assumed errors of the two variables. The uncertainty contained in the MAX-DOAS measurements may be difficult to assess, but, since both MMF and MAPA retrievals are based on the same input data, the associated errors are assumed the same and equal to the mean error provided by MMF and MAPA for each data point.

Table 3

The fraction of the data ( %) that are flagged as valid by MMF and MAPA (individually and combined) for each species.

Species	Inversion algorithm		Combined flagging
	MAPA	MMF
NO${}_{\mathrm{2}}$	29.0	62.4	23.9
HCHO	18.0	82.6	16.8
Aerosols (VIS)	47.6	57.4	33.4
Aerosols (UV)	38.3	54.8	24.8

3.1 Simulated dSCDs

In this section we evaluate the performance of the forward models of MMF and MAPA by intercomparing the simulated trace gas dSCDs of the four species for the entire period. Also, we assess their ability to successfully simulate the slant column densities under different atmospheric (pollution and meteorological) conditions and viewing geometries by comparing the modeled with the measured dSCDs (Fig. ). Each row corresponds to a different trace gas, with the left column presenting the intercomparison results of the modeled dSCDs, while the middle and right columns show the comparison results between the measured dSCDs and the dSCDs simulated by MAPA and MMF, respectively. The data points are colored by the elevation angle, and hotter colors represent dSCDs close to the horizon. A generally better performance of both algorithms is observed for the species retrieved in the VIS range compared to those retrieved in the UV. The modeled slant columns agree well, with Pearson's correlation coefficients and slopes close to unity ($R=\mathrm{0.999}$, $\text{slope}=\mathrm{1.008}$ for NO${}_{\mathrm{2}}$ and $R=\mathrm{0.998}$, $\text{slope}=\mathrm{1.006}$ for O${}_{\mathrm{4}}$ VIS). Additionally, the simulated dSCDs are in good agreement with the measured dSCDs, which is a good indicator for successful profile retrievals. In the case of O${}_{\mathrm{4}}$ (UV), even though the slope and correlation coefficient are similar to O${}_{\mathrm{4}}$ (VIS), a larger scatter is evident, while for HCHO larger deviations from unity in the slopes and correlation coefficients are observed, especially at higher elevation angles. This can probably be explained by the increased noise in the UV spectra compared to the VIS range and also due to the fact that at higher elevation angles the measured differential optical densities are very low, reaching the spectrometer's detection limit. For aerosols in both spectral ranges, discrepancies between the simulated dSCDs of MMF and MAPA may arise due to the variable O${}_{\mathrm{4}}$ scaling factor that is included in MAPA retrievals. This could also be the main driver of the positive bias for low elevation angles that is found in MMF's O${}_{\mathrm{4}}$ dSCDs (especially in the UV), while the results of MAPA are less affected.

Figure 5

Intercomparison of the dSCDs simulated by MMF and MAPA (left column) and comparison of the dSCDs simulated by MAPA (center column) and MMF (right column) against the measured dSCDs. The elevation angles are denoted by different colors (see scale at the bottom).

[Figure omitted. See PDF]

3.2 Averaging kernels

The averaging kernels (AVKs) of a profile retrieval describe the sensitivity of the retrieved state to the true atmospheric state for each altitude layer. The degrees of freedom (DoFs) are mathematically derived as the trace (or sum of the diagonal elements) of the AVK matrix and quantify the number of independent pieces of information gained from the measurements compared to the a priori knowledge . Both the AVKs and the DoFs can be used to characterize the quality of the retrieved profile. Since only OEM-based inversion algorithms are capable of providing AVKs, the results shown here are derived only by MMF. Figure  shows a typical example of the calculated AVKs for each of the retrieved species, including their corresponding DoFs. The median DoFs retrieved by MMF are $\mathrm{3.13}\pm \mathrm{0.32}$ for NO${}_{\mathrm{2}}$, $\mathrm{2.22}\pm \mathrm{0.34}$ for HCHO, $\mathrm{2.73}\pm \mathrm{0.28}$ for aerosols in the VIS, and $\mathrm{2.02}\pm \mathrm{0.32}$ for aerosols in the UV. The averaging kernels illustrate that MAX-DOAS measurements are typically less sensitive for altitudes greater than $\sim \mathrm{2}$ km, as a result of the viewing geometry, and thus, altitudes greater than 3 km are not presented here. That means that the MAX-DOAS measurements under these viewing geometries and with the a priori profiles and covariance matrices used in this study (Sect. ) are adequate for retrieving the extinction and concentration profiles only up to the lowermost $\sim$ 1.5–2 km of the atmosphere with the highest sensitivity closer to the ground. Also, since the photon path increases with wavelength, the MAX-DOAS technique shows higher sensitivity for the species retrieved in the VIS range than in the UV.

Figure 6

A typical example of the retrieved averaging kernels for different altitudes of each species. Hotter colors correspond to altitudes closer to the ground.

[Figure omitted. See PDF]

3.3 Integrated columns

In the past, the trace gas VCDs measured by our MAX-DOAS systems have been derived by dividing the measured dSCDs, only at two elevation angles, 30 and 15${}^{\circ }$, or the mean of the two, with appropriate dAMFs. The dAMFs have been calculated either following the geometrical approximation approach or by deploying RTM simulations taking into account the viewing geometry, the aerosol optical properties, and the instrument's viewing direction relative to the sun . However, in both cases, the actual trace gas profile has not been taken into consideration, introducing, possibly, an additional uncertainty to the measured VCD. This is the first time during the Phaethon's operation that the whole elevation sequence is used in order to retrieve the tropospheric VCDs more accurately. The comparison of the NO${}_{\mathrm{2}}$ and HCHO VCDs that are derived from the integration of the vertical profiles with the VCDs that are calculated using the geometrical approximation can be found in the Appendix .

Figure 7

Time series and scatter plots of the integrated columns for all species retrieved by MMF and MAPA (panel a refers to NO${}_{\mathrm{2}}$, b to HCHO, c to aerosols in the visible range, and d to aerosols in the UV). The parameters of the orthogonal distance regression, i.e., slope ($S$), offset ($O$), number of points ($N$), and Pearson's correlation coefficient ($R$), are shown in different colors for each azimuth viewing direction. The text in black color represents the consolidated statistics for all azimuth directions. The dashed black line represents the $\mathrm{1}:\mathrm{1}$ line.

[Figure omitted. See PDF]

In Fig.  the time series of the integrated columns of all retrieved species (i.e., AODs for aerosols and VCDs for trace gases) are presented, as well as comparisons between MMF and MAPA. The statistics of the comparisons, i.e., slope ($S$), offset ($O$), number of points ($N$), and Pearson's correlation coefficient ($R$), are shown in different colors for each azimuth viewing direction. The text in black color represents the consolidated statistics for all azimuth directions. No clear azimuth dependence of the retrieved columns is observed for the trace gases. However, for aerosols, especially in the UV, significant differences in the regression slopes appear for the different azimuths. It should be noted that in these comparisons a variable O${}_{\mathrm{4}}$ scaling factor has been used for the MAPA retrievals, and since no scaling factor has been applied to the MMF retrievals, differences in the AODs between the two algorithms are expected. The number of elevation sequences at 220${}^{\circ }$ azimuth is always larger compared to the other azimuth directions, because the instrument was configured to record spectra only at this particular direction for approximately 1 month in the beginning of its operation. The comparison shows that the NO${}_{\mathrm{2}}$ VCDs derived by MMF and MAPA are in very good agreement, with slopes and correlation coefficients close to unity (ranges: $\mathrm{0.999}\le S\le \mathrm{1.024}$ and $\mathrm{0.977}\le R\le \mathrm{0.985}$). Similar results were obtained for all azimuth directions with $S=\mathrm{1.009}$ and $R=\mathrm{0.982}$. In the case of HCHO, despite the good correlation ($R=\mathrm{0.927}$), notable deviations from unity in the slope are observed for all azimuth directions. MAPA systematically reports larger VCDs than MMF for higher HCHO concentrations, while the opposite behavior is observed for low HCHO loads, indicating that further investigation is required. This behavior could be explained by the increased spectral noise in the UV that leads to discrepancies between the HCHO dSCDs simulated by the forward models of MMF and MAPA (see discussion in Sect. ) and due to an invariant a priori profile during the year. The sensitivity of the MAX-DOAS decreases with altitude and it is very limited at altitudes above $\sim \mathrm{2.5}$ km for the species measured in the VIS spectral range or even lower ($\sim \mathrm{1.5}$ km) for the species in the UV (Fig. ). For NO${}_{\mathrm{2}}$, this is generally not a problem since the total column is dominated by the concentration in the lower layers of the troposphere (see also discussion in Sect. ). However, HCHO can be vertically extended at higher altitudes, where the sensitivity of the MAX-DOAS is low. In the case of HCHO, MMF is more prone to result in the a priori profile, while MAPA retrievals become more unstable. Thus, the vertical profiles of MAPA are expected to have greater variability. Concerning aerosols, the comparison of the retrieved AODs reveals better agreement at 477 nm ($R=\mathrm{0.856}$) than at 360 nm ($R=\mathrm{0.700}$), with larger scatter and more outliers compared to the trace gas VCDs. As already mentioned, this is mainly attributed to the O${}_{\mathrm{4}}$ scaling factors that are used in MAPA retrievals. More details about the effect of the O${}_{\mathrm{4}}$ scaling factor on the retrieved AODs and the trace gas VCDs can be found in Appendix .

3.4 Surface concentrations

The surface concentration is defined as the trace gas amount at ground level. However, the profile parameterization used within MAPA allows for the retrieval of lifted trace gas layers for a shape parameter greater than 1, which leads to a value of zero for the concentration at the surface. For these cases, the comparison with in situ measurements or surface concentrations retrieved by an OEM-based algorithm will be low-biased. Thus, in the following sections, the term “surface concentration” will refer to the mean “near-surface concentration”, i.e., the average concentration below 200 m for both MMF and MAPA, rather than the concentration directly at the ground. Figure  shows the time series of the near-surface NO${}_{\mathrm{2}}$ and HCHO concentrations derived by MMF and MAPA and the corresponding scatter plots. The comparisons of the surface values are similar to the comparisons of the tropospheric VCDs (shown in Fig. ) with $\text{slope}=\mathrm{1.118}$, $R=\mathrm{0.919}$ for NO${}_{\mathrm{2}}$ and $\text{slope}=\mathrm{1.295}$, $R=\mathrm{0.855}$ for HCHO, but more outliers are present. In the case of HCHO, the surface concentrations derived for the 142${}^{\circ }$ azimuth direction show larger differences compared to the other directions, while this is not clear for NO${}_{\mathrm{2}}$. These discrepancies are possibly related to the fact that for this azimuth, the elevation angle of 1${}^{\circ }$ was not included in the analysis (see Sect. ), which may have influenced the retrieved surface concentrations.

Figure 8

Time series and scatter plots of the near-surface concentrations of NO${}_{\mathrm{2}}$ (a) and HCHO (b) derived by MMF and MAPA. The parameters of the orthogonal distance regressions are presented as in Fig. .

[Figure omitted. See PDF]

3.5 Seasonal mean vertical profiles

Figure  shows the seasonal mean NO${}_{\mathrm{2}}$, HCHO, and aerosol extinction vertical profiles at Thessaloniki retrieved by MMF (cyan) and MAPA (magenta) during the 12 months considered in this study. Each row represents the vertical profiles of a specific species, and each column corresponds to a different season. The shaded areas represent the standard deviation around the mean for each layer, and they illustrate the seasonal variability of the vertical profiles. For NO${}_{\mathrm{2}}$ both algorithms report profiles that are decreasing with altitude for all seasons. Compared to MAPA, MMF reports slightly lower NO${}_{\mathrm{2}}$ concentrations below 1 km (yet within the range of variability) and slightly higher above 1 km. For all seasons, the variability of MMF's profiles between 1 and 2 km is larger, probably due to the increased contribution of the a priori profile under certain conditions (e.g., high aerosol load or fog close to the ground), where the sensitivity of the MAX-DOAS is lower. However, the seasonal mean profiles of both algorithms indicate that most of the NO${}_{\mathrm{2}}$ content lies within the first $\sim \mathrm{500}$ m. NO${}_{\mathrm{2}}$ originates mainly from direct, local emissions close to the ground (e.g., road transport emissions). Additionally, its lifetime in the PBL is short, typically a few hours depending on the season e.g.,, and as a result, higher NO${}_{\mathrm{2}}$ amounts are expected at lower altitudes in the troposphere. Both algorithms retrieve higher concentrations close to the surface during the cold period, which can be mainly attributed to enhanced NO${}_{\mathrm{2}}$ emissions near the ground (e.g., from domestic heating sources) and to reduced photolysis rates due to weaker solar radiation.

An opposite seasonal variation is observed for HCHO, with higher concentrations reported by both algorithms during summer (consistent with the VCDs, shown in Fig. ). The profile shapes of MMF and MAPA agree reasonably well. In summer, the larger retrieved concentrations are probably due to the increased emissions of VOCs, whose oxidation produces HCHO. According to , biogenic emissions are expected to peak during summer, while the anthropogenic emissions do not show a clear seasonal variation in Thessaloniki. Therefore, the observed HCHO seasonality is mainly attributed to the enhanced biogenic emissions from vegetation in summer. VOCs are generally well mixed and have longer life times , hence, larger HCHO amounts are expected at higher altitudes during the warm season. MMF's profiles peak at a slightly higher altitude ($\sim \mathrm{800}$ m) than MAPA's ($\sim \mathrm{500}$ m) and decrease with a slightly higher rate and less variability for altitudes above 1 km. However, such differences, especially at higher altitudes, are to some extent expected, since the sensitivity of the MAX-DOAS decreases rapidly with altitude for the species that are measured in the UV (Fig. ). This means that concentrations at high altitudes are strongly constrained by the a priori profile in the retrievals of MMF. Also, parameterized algorithms (such as MAPA) have the tendency of becoming unstable when the sensitivity is low .

For aerosols, the largest differences in the vertical profiles of MMF and MAPA are found in the VIS range and especially in summer and autumn. MMF yields more structured aerosol extinction profiles for altitudes between 1 and 2 km, while MAPA reports smoother, exponentially decreasing profiles. Such differences are not found in the UV retrievals during summer and autumn. It should be noted that larger discrepancies among different inversion algorithms for the species retrieved in the VIS compared to the UV have also been reported in other studies e.g.,. At higher altitudes, both algorithms report greater aerosol concentrations in summer than the winter, in both the UV and VIS. Similar results were found for Thessaloniki by using seasonal mean vertical profiles measured by a lidar system. In contrast, near the surface, aerosol concentrations are highest in winter and lowest in summer. This pattern can be mainly attributed to the shallower PBL in Thessaloniki during winter and autumn that shrinks to $\sim \mathrm{1}$ km due to the weaker solar radiation and lower air temperatures. In the UV during winter, MAPA retrieves larger aerosol concentrations close to the ground, which decrease more rapidly with altitude than MMF. However, as already discussed for HCHO, the sensitivity of the MAX-DOAS in the UV is very limited at higher altitudes, and the profiles of MMF are driven towards the a priori profile. Another contributing factor for the differences in the aerosol extinction profiles between the two algorithms might be the variable O${}_{\mathrm{4}}$ scaling factor that is used in MAPA retrievals, while no scaling factor is applied in MMF. The effect of the O${}_{\mathrm{4}}$ scaling factor on the AOD is presented in the Appendix .

Figure 9

The seasonal mean vertical profiles of NO${}_{\mathrm{2}}$, HCHO, and aerosol extinction in the VIS and UV retrieved by MMF and MAPA. Each row (1–4) represents the vertical profiles of different species along the four seasons (columns 1–4). The shaded areas represent the standard deviation around the mean for each layer, and they illustrate the seasonal variability of the profiles.

[Figure omitted. See PDF]

4 Validation

In this section we present the validation results of the products retrieved by the MAX-DOAS profile analysis against ancillary data measured by other reference co-located instruments. Vertical profiles of the aerosol extinction measured by a co-located lidar system are used to validate the aerosol vertical profiles retrieved by the MAX-DOAS, while the AODs in the UV and VIS range are compared with those measured by a sun photometer and a spectrophotometer. The NO${}_{\mathrm{2}}$ near-surface concentrations are compared with in situ surface measurements, but since no other sources of HCHO data are available, the MAX-DOAS-derived vertical profiles, columns, or surface concentrations cannot be validated.

4.1 Aerosol extinction profiles

The AOD values at 477 and 360 nm retrieved by the MAX-DOAS are compared with the AOD measured by the co-located CIMEL sun photometer and the Brewer spectrophotometer. Quasi-simultaneous (within $\pm \mathrm{15}$ min) measurements were found, and the AODs at 477 and 360 nm were calculated using the Ångström exponent between 380 and 500 nm and the AOD at these wavelengths derived by the CIMEL sun photometer. Since the Brewer spectrophotometer's wavelength range spans up to 365 nm, climatological monthly mean Ångström exponent values, calculated from the CIMEL data, have been used to extrapolate the AOD to 477 nm. Figure  shows the time series of all AOD data at 477 and 360 nm (not just the quasi-simultaneous) retrieved by the three systems. The CIMEL sun photometer was not operating for approximately 4 months during the summer of 2020 due to a delay in its scheduled annual maintenance and calibration. AOD data derived by the Brewer spectrophotometer are available until January 2021.

Figure 10

Time series of all available AOD data at 477 and 360 nm retrieved by the MAX-DOAS system, the Brewer spectrophotometer, and the CIMEL sun photometer.

[Figure omitted. See PDF]

Since MMF and MAPA rely on their own individual flagging schemes in order to ensure that the retrieved products are of high quality, we investigate the effect of applying different flagging schemes to the data, which are listed in Table .

Table 4

The flagging schemes that are applied to the retrieved products.

Flagging scheme	Description
No. 1	Data are flagged as valid by the flagging algorithm of MMF
No. 2	Same as scheme no. 1, but for the flagging algorithm of MAPA
No. 3	Data that are flagged as warning by either MMF or MAPA are also considered valid
No. 4	Data are flagged as valid by both MMF and MAPA

Scheme nos. 1 and 2 correspond to the default own flagging algorithms of MMF and MAPA, scheme no. 4 is expected to provide data of maximum quality since data are designated as valid by both algorithms, while scheme no. 3 rejects the error-flagged data but treats the warnings raised by MMF or MAPA as valid data. Figure  shows the comparison between the common AOD data derived by the CIMEL sun photometer, Brewer spectrophotometer, and the MAX-DOAS at 360 and 477 nm. Each column of the figure corresponds to a different flagging scheme as described in Table . Figure  graphically presents the statistics of the linear regressions (i.e., slope, offset, number of points, and Pearson's correlation coefficient) between the reference instruments and the MAX-DOAS. The panels (a)–(d) correspond to different flagging schemes, as Fig. .

Figure 11

Scatter plots of the AOD retrieved by the MAX-DOAS data analyzed by MMF and MAPA against the CIMEL sun photometer (a at 477 nm, c at 360 nm) and the Brewer spectrophotometer (b at 477 nm, d at 360 nm). Each column represents data that are flagged as valid according to the flagging schemes of Table .

[Figure omitted. See PDF]

Figure 12

Graphical representation of the linear regression parameters (slope, offset, correlation coefficient, and number of data) of the comparison between the AOD derived from MAPA and MMF against the AOD from the CIMEL sun photometer (a, c) and the Brewer spectrophotometer (b, d) at 477 nm (a, b) and 360 nm (c, d) for each flagging scheme (nos. 1 to 4).

[Figure omitted. See PDF]

The comparison results of the MAX-DOAS against the CIMEL sun photometer are slightly different than those with the Brewer spectrophotometer. This can probably be explained by the fact that only a few collocated measurements are available and in different periods for the two reference instruments (Fig. ). In the case of the AOD at 477 nm most of the outliers are filtered out when the flagging scheme no. 4 is applied, and the best agreement is observed between the reference instruments and the MAX-DOAS, for both MMF and MAPA, with similar correlation coefficients (0.79 for the CIMEL sun photometer and 0.81 for the Brewer spectrophotometer). Compared to the CIMEL sun photometer, MAPA seems to perform slightly better than MMF when each algorithm considers its own flagging, with correlation coefficients of 0.70 and 0.50, respectively (MAPA for scheme no. 2 and MMF for scheme no. 1). However, compared to the Brewer spectrophotometer, both algorithms show very similar correlation coefficients (0.78 and 0.77). Swapping the flags between MMF and MAPA leads to worse agreement and more outliers. The results of scheme no. 3 indicate that some of the warning-flagged data are of lower quality and should be treated with caution.

In the case of the AOD at 360 nm, the effect of the flagging schemes is different. Here, most of the outliers are eliminated, and the best overall agreement is achieved for scheme no. 2 (i.e., when MAPA's individual flagging algorithm is applied). This behavior is observed for both MMF and MAPA, indicating that the flagging algorithm of MAPA performs better than that of MMF in the UV. The correlation coefficients with the CIMEL sun photometer data are 0.72 for MAPA and 0.70 for MMF, and with the Brewer spectrophotometer data they are 0.72 for MAPA and 0.78 for MMF. The flagging scheme no. 4 removes even more data (as expected) but does not improve the comparisons. The effect of the warning-flagged data (scheme no. 3) is more apparent in the case of the UV, and the results suggest that they should not be considered valid. The AOD derived from the MAX-DOAS, in both the UV and the VIS range, is, generally, underestimated compared to the AOD measured by the CIMEL sun photometer and the Brewer spectrophotometer. This finding is consistent with other studies e.g.,. However, it should be noted that the AODs derived by the MAX-DOAS and the reference instruments may not always refer to the same physical quantity. The CIMEL sun photometer and the Brewer spectrophotometer use direct-sun observations to retrieve the total column amount of the aerosol extinction, while the MAX-DOAS sensitivity decreases rapidly with altitude (Sect. ), and the derived AOD corresponds mainly to the lowermost tropospheric aerosol (partial AOD). Additionally, the vertical profiles retrieved by OEM-based algorithms are biased towards the a priori profile at higher altitudes (Fig. ), leading to deviations from the true profile, meaning that aerosol layers above $\approx \mathrm{2}$ km cannot be reliably retrieved . Discrepancies in the AODs between the instruments are expected, usually when aerosols are present at altitudes greater than 2 km, contributing to the total AOD, but not detected by the MAX-DOAS.

During the whole period of this study only a few lidar measurements of the aerosol extinction profile are available, so the true state of the aerosol profile is generally not known. Additionally, synchronous measurements between the lidar and the MAX-DOAS are even fewer, so an in-depth validation of the aerosol profiles retrieved by the MAX-DOAS needs further investigation. In this section we present the comparison of four profiles retrieved by the two systems within $\pm \mathrm{30}$ min (using the MAX-DOAS profile that is closest in time with the lidar measurement) and which are indicative for the period of study. An important issue that arises in the validation of the MAX-DOAS vertical profiles is that usually the validator (in this case the lidar system) allows the detection of aerosol layers in a much higher vertical resolution than the MAX-DOAS. When the true aerosol profile state is actually known and in order to compare meaningful results, the lidar profiles need to be smoothed (i.e., degraded to the sensitivity of the MAX-DOAS). Only the OEM-based algorithm provides information that can be used to smooth the lidar profiles. The information about the sensitivity is quantified by the averaging kernel according to : 2$${\mathit{x}}_{\mathrm{s}}={\mathit{x}}_{\mathrm{a}}+\mathbf{A}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right),$$ where ${\mathit{x}}_{\mathrm{s}}$ is the smoothed lidar profile, ${\mathit{x}}_{\mathrm{a}}$ and $\mathbf{A}$ are the a priori profile and the averaging kernel of the OEM-based retrieval, and $\mathit{x}$ is the initial lidar profile. Deviations of the smoothed lidar profile at each altitude depend on the a priori profile and the sensitivity of the MAX-DOAS at this altitude . Since the sensitivity of the MAX-DOAS decreases with altitude, the application of the averaging kernels is expected to smooth the true profiles towards lower altitudes. However, since MAPA does not quantify the sensitivity, a similar smoothing cannot be performed for MAPA retrievals, and thus, the aerosol extinction profiles are directly compared with the initial lidar profiles. Another point that should be noted is the differences in the operational principles of the two instruments. The lidar retrieves the vertical profile from the air mass that is located overhead, while the MAX-DOAS scans through different air masses along the line of sight of the telescope during an elevation sequence . Its effective horizontal distance is of the order of a few kilometers and increases at elevation angles close to the horizon. Thus, differences in the retrieved extinction profiles are expected, especially at locations with large horizontal inhomogeneities of aerosols. As already mentioned (Sect. ) a constant climatological lidar ratio of 50 sr was assumed for the channel of 532 nm and was applied to the backscatter profiles in order to retrieve the extinction, which may also result in uncertainties of the validator's product. So, in this study the comparisons are focused on the shape of the profiles and the retrieved aerosol layer heights rather than on the absolute values of the aerosol extinction.

Figure 13

Four cases of aerosol vertical profiles measured by the MAX-DOAS (cyan for MMF and magenta for MAPA) and the lidar (black for the original and yellow for the smoothed profiles). The shaded area represents altitudes where the lidar is not capable of retrieving the aerosol extinction profile accurately due to the overlap effect.

[Figure omitted. See PDF]

Figure  presents the comparison of four aerosol extinction profiles in the VIS retrieved by the MAX-DOAS and the lidar. The lidar profile is trustworthy only above a certain altitude (approximately 0.6 km) owing to the geometry of the telescope and the emitted laser beam, which prevents a fraction of the backscattered radiation from reaching the detector at altitudes close to the surface (overlap effect). Thus, the aerosol extinction retrieved by the lidar below 0.6 km is not presented. Since the MAX-DOAS profile retrievals in the UV are sensitive at altitudes closer to the ground (see discussion in Sect. ), where the lidar system is not, the profiles for 360 nm are excluded from the analysis. Less information content obtained for species measured in the UV was also apparent in other studies e.g.,. In Fig.  the yellow line corresponds to the smoothed lidar profile, which is degraded to the sensitivity of the MAX-DOAS according to Eq. (). In general, the aerosol vertical profiles retrieved by the MAX-DOAS can realistically estimate the shape of the true state, even though some differences appear between MMF and MAPA. The agreement between the shape of the MMF profiles and the smoothed lidar profiles is much better than for MAPA. This is expected since the initial lidar profile is degraded to the sensitivity of the MAX-DOAS using the averaging kernels derived by MMF and the a priori profile of the retrieval. On 4 and 5 June the aerosol load is low, and both algorithms report similar profiles. In both cases the MAX-DOAS profiles successfully fit the shape of the true profile. However, on 5 June the shape of the profiles is similar only up to $\approx \mathrm{2}$ km. Aerosol layers between 2 and 4 km that are detected by the lidar cannot be well retrieved by the MAX-DOAS, due to its very limited sensitivity at these altitudes (Fig. ). In the other two cases MMF and MAPA provide different profiles. Discrepancies in the aerosol extinction profiles retrieved in the VIS by MMF and MAPA are also found in the seasonal mean vertical profiles (Sect. ), especially during summer. On 21 July MMF reports a more structured aerosol profile with a distinct aerosol layer at about 1.3 km, while MAPA reports a smoother profile with a thick layer spanning from the surface up to $\approx \mathrm{2}$ km. On 28 August, even though MMF and MAPA report profiles of different shapes below $\approx \mathrm{1}$ km, the profiles agree reasonably well with the lidar profile. At higher altitudes MMF is biased towards the a priori profile due to the limited sensitivity, while the aerosol extinction retrieved by MAPA decreases rapidly. Despite the observed differences, the results of the comparisons are promising, indicating that the analysis of the MAX-DOAS data can provide a generally good estimation of the vertical aerosol extinction profiles over Thessaloniki. However, further investigation is required in order to assess the differences in the aerosol profiles provided by the two systems, but also between the two inversion algorithms, when more collocated measurements become available. Such studies will be further facilitated with a new lidar system with improved overlap height that is currently under development, which will allow the retrieval of the aerosol profiles at altitudes closer to the ground, where the MAX-DOAS shows higher sensitivity.

4.2

NO${}_{\mathrm{2}}$ surface concentration

In Fig.  we present a comparison of near-surface concentrations derived from the MAX-DOAS data (i.e., the average concentration 200 m) with in situ NO${}_{\mathrm{2}}$ measurements. The small dots represent the hourly mean values, while the solid lines refer to the daily mean concentrations. The comparison is only performed for the hourly mean concentrations derived by the two systems, while the daily mean concentrations are shown only for a qualitative comparison. The MAX-DOAS hourly mean values are horizontally averaged from all azimuth viewing directions. A dataset from June 2020 to March 2021 (about 10 months) is considered in this study, in which both MAX-DOAS and in situ measurements were available. The MAX-DOAS reports systematically lower NO${}_{\mathrm{2}}$ concentrations than the in situ by $\sim$ 55 %–60 %. Since the concentrations retrieved by the MAX-DOAS are averaged along a horizontal path of a few kilometers, which may extend over the bay, whereas the in situ data refer to a specific location (point measurements), the MAX-DOAS is generally expected to report lower values from the air quality station, which is also occasionally affected by local emissions. Differences in the retrieved concentrations may also arise due to the slightly different altitudes of the measurement sites. Similar results have been found in other studies e.g.,. Both MMF and MAPA perform well for the retrieval of the NO${}_{\mathrm{2}}$ near-surface concentrations. Even though the in situ site is located opposite to the MAX-DOAS system's azimuth viewing directions (see Fig. ), the correlation coefficients of both algorithms are still high, suggesting that strong NO${}_{\mathrm{2}}$ horizontal inhomogeneities are less likely to occur in Thessaloniki, at least during the period of study. The effect of the different flagging schemes is not as strong as for aerosols, except for MAPA when MMF's flagging (scheme no. 1) is applied to the data. The performance of MMF is slightly better than MAPA's with fewer outliers and higher correlation coefficients, even though it reports a much larger fraction of data as valid (Table ). The results of flagging scheme no. 3 indicate that the warning-flagged data could also be considered valid. This could be explained by the fact that a large part of the flagged data are related to the effects of clouds. As shown in previous studies e.g.,, under most cloud conditions (except for fog and optically thick clouds), the trace gas vertical profiles, columns, and near-surface concentrations can still be retrieved, while the aerosol retrievals are more strongly affected.

Figure 14

(a) Time series of hourly mean NO${}_{\mathrm{2}}$ surface concentrations derived from MAX-DOAS by MMF (cyan) and MAPA (magenta) and from in situ measurements (black). The solid lines correspond to the time series of daily means. Panels (b) and (c) are the corresponding scatter plots colored by the four flagging schemes that are applied to the MAX-DOAS data, and panel (d) shows the statistics of the comparisons.

[Figure omitted. See PDF]

5 Summary and conclusions

In this study we have retrieved vertical profiles of aerosols, NO${}_{\mathrm{2}}$, and HCHO for the first time in Thessaloniki, Greece, using MAX-DOAS observations by applying an OEM-based inversion algorithm (MMF) and a parameterized algorithm (MAPA). Their performance is evaluated by intercomparing the dSCDs simulated by the forward models, the integrated columns (i.e., VCDs for trace gases and AODs for aerosols), the trace gas near-surface concentrations, and the seasonal mean vertical profiles derived by the two algorithms. The products that are retrieved by the inversion analysis of MAX-DOAS measurements using MMF and MAPA are compared with ancillary data measured by other reference instruments. This study provides the basis for future research activities, e.g., the investigation of the spatio-temporal variability of trace gas and aerosol profiles over Thessaloniki.

The tropospheric column densities of NO${}_{\mathrm{2}}$ are in excellent agreement (slope very close to unity and $R=\mathrm{0.982}$), while for HCHO, even though a generally good correlation is found ($R=\mathrm{0.927}$), deviations from unity in the slopes are observed, which can be attributed to discrepancies between the HCHO dSCDs simulated by the forward models of MMF and MAPA and the limited sensitivity of the MAX-DOAS in the UV, especially at higher altitudes. Concerning aerosols, a better agreement between MMF and MAPA is found for the AOD at 477 nm than at 360 nm due to the increased SNR in the VIS range and the stronger effect of the O${}_{\mathrm{4}}$ scaling factor on the retrieved AODs in the UV. No clear azimuth dependence is observed for any of the retrieved species. The seasonal mean vertical profiles retrieved by MMF and MAPA are generally in good agreement. The largest discrepancies are found for the aerosol extinction profiles in the VIS and especially during summer.

The AODs retrieved by the MAX-DOAS are validated by comparison with measurements of a CIMEL sun photometer and a Brewer spectrophotometer. Four flagging schemes were applied to the MAX-DOAS-derived data, and their effect on different products is evaluated. However, no robust conclusion could be drawn about which flagging algorithm shows an overall better performance. A generally good qualitative agreement is found for both VIS and UV wavelengths (with correlation coefficients up to 0.8). The negative bias that is observed from the reference instruments is probably mostly due to the limited sensitivity of the MAX-DOAS in retrieving aerosol information at higher altitudes, especially in the UV. The results also indicate that using an intersected dataset derived by applying both flagging algorithms to the data improves the agreement of AODs at 477 nm; however, a similarly good agreement is not observed in the UV, where the flagging algorithm of MAPA performs better. Four cases of aerosol extinction vertical profiles at 477 nm are compared with profiles measured by a co-located lidar system. The MAX-DOAS was found to provide a generally good estimation of the shape of the profile. The NO${}_{\mathrm{2}}$ near-surface concentrations are compared with in situ NO${}_{\mathrm{2}}$ observations, where the effect of the different flagging schemes is not found to have such a strong impact as for aerosols. The concentrations from both MMF and MAPA are in good agreement with the in situ measurements in terms of variability but are highly biased by approximately 60 %. MMF shows a slightly better performance ($R=\mathrm{0.78}$) compared to MAPA ($R=\mathrm{0.73}$).

The effect of the O${}_{\mathrm{4}}$ scaling factor is also investigated by comparing the integrated columns of MMF and MAPA and also by comparing the AODs derived by MAPA for different values of the scaling factor with AODs measured by the CIMEL sun photometer and the Brewer spectrophotometer. The effect of the O${}_{\mathrm{4}}$ scaling factor has a stronger impact on aerosols than on trace gases (where the effect is minor). The fixed value of 0.8 for the scaling factor, which is supported by many studies, does not seem to be suitable for the measurements at Thessaloniki.

Appendix A

Effect of the O${}_{\mathrm{4}}$ scaling factor

The O${}_{\mathrm{4}}$ scaling factor (SF) was introduced by in order to remove the systematic discrepancies appearing between measured and simulated O${}_{\mathrm{4}}$ dSCDs. Uncertainties of the O${}_{\mathrm{4}}$ cross sections and/or its temperature and pressure dependence, aerosol optical properties, and RTM errors have been suggested as possible causes for these discrepancies. Several studies have confirmed the idea of the O${}_{\mathrm{4}}$ scaling factor and have shown that applying a SF (commonly using a value between 0.75 and 0.9) is indeed necessary . However, other studies have not supported this requirement . Although the need for an O${}_{\mathrm{4}}$ scaling factor for retrieving aerosol information from MAX-DOAS measurements has been extensively discussed , its physical mechanism is not understood and still remains an unresolved issue.

Figure A1

Scatter plots of the integrated columns of NO${}_{\mathrm{2}}$, HCHO, and AOD at 477 and 360 nm (a–d, respectively) retrieved by MMF and MAPA for various O${}_{\mathrm{4}}$ scaling factors (columns 1–4).

[Figure omitted. See PDF]

Since MAPA provides the option of scaling the modeled O${}_{\mathrm{4}}$ dSCDs, we have investigated the effect of the SF on the comparisons between the products of the two profiling algorithms and between AOD derived by MAPA and the reference instruments. We selected three fixed values (i.e., 1, 0.9, and 0.8, referred to hereafter as SF1.0, SF0.9, and SF0.8) and a variable SF (SFvar). Figure  shows the effect of the O${}_{\mathrm{4}}$ SF on the comparison of trace gas VCDs and AOD derived by MAPA and MMF (same as Fig.  without accounting for the different azimuth directions). Like in Fig. , the regression results are based on an ODR, and also the retrievals of aerosols in the UV are based on the flagging algorithm of MAPA. Since MMF does not take into account a scaling for O${}_{\mathrm{4}}$, the use of different SFs in MAPA leads to substantial differences in the regression slopes and correlation coefficients when comparing the AODs of the two algorithms, for both 360 and 477 nm. The closest to unity slopes and the highest correlation coefficients are found, as expected, when no scaling factor is applied (SF1.0) for both wavelengths. The slopes of the fitting are 1.18 and 1.07 for the AODs at 477 and 360 nm, respectively, while the correlation coefficients are $\sim \mathrm{0.89}$ and $\sim \mathrm{0.88}$. Especially for the AOD in the UV, the agreement between MMF and MAPA for the variable SF substantially declines and the scatter increases. The worst results appear for SF0.8, also leading to substantial reduction in the reported valid data. The use of a scaling factor does not seem to affect the retrieved VCDs for NO${}_{\mathrm{2}}$ and HCHO, at least as it concerns the slope and correlation coefficient of the regression, but there is some effect on the number of the data reported as valid, especially for SF0.8. Opposite to aerosols, the best correlation of the retrieved NO${}_{\mathrm{2}}$ and HCHO columns is achieved when using the SFvar instead of the SF1.0.

Figure A2

Comparison of the AODs at 477 nm (a, b) and 360 nm (c, d) derived by MAPA with the AODs measured by the CIMEL sun photometer (a, c) and the Brewer spectrophotometer (b, d) for different O${}_{\mathrm{4}}$ SFs (columns 1–4).

[Figure omitted. See PDF]

Figure  presents the comparison of the AODs at 360 and 477 nm retrieved by the MAX-DOAS (using MAPA) with the AODs calculated by the CIMEL sun photometer and the Brewer spectrophotometer for different O${}_{\mathrm{4}}$ SFs. For consistency, the individual flagging algorithm of MAPA is used for the retrievals in both the UV and VIS. The differences in the slopes and correlation coefficients among the different SFs are rather small, with the former dominated mainly by the noise in the measurements, as discussed in Sect. . The results of SF1.0 (i.e., no scaling factor) show a similar performance, with the SFvar indicating that the most suitable O${}_{\mathrm{4}}$ SF value for Thessaloniki would be closer to unity. As for the AOD, the SF0.8 flags a larger fraction of the data as invalid compared to the other SFs in all cases, and the agreement between the MAX-DOAS and the reference instruments gets worse with decreased correlation coefficients and larger offsets and slopes. Hence, the SF0.8 that is supported by many studies for achieving better agreement between the MAX-DOAS and sun photometers is found to be too small for the profiling of the MAX-DOAS measurements at Thessaloniki.

This is also supported by the frequency distribution of the fitted O${}_{\mathrm{4}}$ SFs in the UV and VIS ranges, shown in Fig. a. The median O${}_{\mathrm{4}}$ SF fitted by MAPA is $\mathrm{0.964}\pm \mathrm{0.125}$ and $\mathrm{0.964}\pm \mathrm{0.115}$ at 477 and 360 nm, respectively. The histograms indicate that for most elevation sequences a scaling factor close to unity is required to bring measured and simulated O${}_{\mathrm{4}}$ dSCDs into agreement. The SF0.8 seems to be too low for the retrievals in Thessaloniki for both spectral ranges. However, it should be noted that an apparent seasonal pattern in the fitted O${}_{\mathrm{4}}$ SF is observed at both 477 and 360 nm, shown in Fig. b. In order to remove any possible effects of the seasonal variability of the SZA, only the O${}_{\mathrm{4}}$ SFs for which $\mathrm{65}{}^{\circ }<\text{SZA}<\mathrm{75}{}^{\circ }$ are presented. The maximum O${}_{\mathrm{4}}$ SF values are reported in August ($\sim \mathrm{1.02}$) at 477 nm and in September ($\sim \mathrm{1.05}$) at 360 nm, while the minimum O${}_{\mathrm{4}}$ SF values are found in February–March for both wavelengths (i.e., $\sim \mathrm{0.86}$ at 477 nm and $\sim \mathrm{0.91}$ at 360 nm). This seasonal variability could partly be explained by the temperature dependence of the O${}_{\mathrm{4}}$ absorption. Since the absorption is stronger at lower temperatures, higher O${}_{\mathrm{4}}$ SFs are generally expected during summer than in winter. The seasonal pattern could also be related to the similar seasonal variability of the AOD. In general, higher AODs are observed over Thessaloniki in summer than in winter . The O${}_{\mathrm{4}}$ SFs for the two wavelengths show a similar but not identical seasonality, and thus, further investigation is required when more MAX-DOAS data become available.

Figure A3

Frequency distribution of the fitted O${}_{\mathrm{4}}$ scaling factor at 477 and 360 nm (a) and its seasonal variability (b). The error bars represent the standard deviation around the monthly averages.

[Figure omitted. See PDF]

Appendix B Comparison with the VCDs calculated using the geometrical approximation

In this section, we compare the trace gas VCDs obtained from the integration of the vertical profiles retrieved by MMF and MAPA with those derived by other methods used in the past at Thessaloniki, in order to establish a link with the VCDs reported in former studies e.g.,. As already discussed in Sect.  the main product that is derived from DOAS spectral analysis is the trace gas dSCD at different viewing directions. However, the conversion of dSCDs to VCDs is usually challenging. The easiest approach that has been adopted for several years is the so-called geometrical approximation . This methodology considers only the geometric light path through the troposphere for the attenuation of radiance at an elevation angle α, and thus the AMF that is required for the calculation of the VCD can be geometrically approximated by B1$${\text{AMF}}_{\mathit{\alpha }}={\displaystyle \frac{\mathrm{1}}{\sin \left(\mathit{\alpha }\right)}}.$$

This approach has been proven potentially appropriate when higher elevation angles are used (typically 30${}^{\circ }$ and/or 15${}^{\circ }$) under low-aerosol conditions . However, for a more accurate calculation of the true AMF, several other parameters must be taken into account, such as the solar position, viewing geometry, ground albedo, wavelength, and aerosol properties . High AOD values are not infrequent in Thessaloniki, especially during summer e.g.,. Thus, for a more accurate calculation of the VCD, calculated dAMF LUTs based on RTM simulations using several parameters (such as the SZA, the AOD, the elevation angle, and the azimuth angle relative to the solar azimuth), appropriate for the fitting windows of the spectrometers used in that study. Since in the current study a different spectral range was used (Sect. ), these LUTs could not be used here.

Tropospheric dSCDs measured at elevation angles of 30 and 15${}^{\circ }$, relative to the zenith, are converted to VCDs by applying the geometrical approximation, and their average is calculated and used when they agree to at least within 50 %. This filtering is necessary since these dSCDs (especially at 30${}^{\circ }$) are much smaller than those measured at lower elevation angles and are associated with larger fitting errors (especially for HCHO). For the VCDs derived from the profiles the default flagging scheme was applied, i.e., when they are flagged as valid by both MMF and MAPA (scheme no. 4 of Table ). In this comparison, no discrimination of the sky conditions (aerosol and cloud) is made.

The comparison for NO${}_{\mathrm{2}}$ and HCHO is shown in Fig. . The number of available data for HCHO is much lower than for NO${}_{\mathrm{2}}$ due to the fewer valid profiles retrieved by the inversion algorithms (Table ) and also because HCHO dSCDs derived from this MAX-DOAS system contain more noise, due to weak SNR in the UV, leading to larger differences between the VCDs obtained from 30 and 15${}^{\circ }$ elevation angles. For NO${}_{\mathrm{2}}$, the profile-derived VCDs compare well to the VCDs obtained from the geometrical approximation with high correlation coefficients ($R=\mathrm{0.96}$ and $R=\mathrm{0.93}$ for MMF and MAPA, respectively). For HCHO the correlation coefficients are lower ($R=\mathrm{0.93}$ for MMF and $R=\mathrm{0.89}$ for MAPA). The discrepancies in the results between MMF and MAPA are consistent with the results of Fig.  and are discussed in Sect. . NO${}_{\mathrm{2}}$ is typically located at lower altitudes than HCHO, mainly because it is produced by emission sources close to the surface (see discussion in Sect. ). Other studies have shown that the error of the geometric VCDs is usually less than 20 % compared to the integrated vertical profiles for trace gases that are located below 1 km, while for trace gases that are located at higher altitudes the geometrical approximation is less accurate since the effect of aerosols becomes more important e.g.,. This finding is consistent with the results of Fig. , where lower correlation coefficients are found for HCHO as HCHO layers can be vertically extended to higher altitudes, making the geometrical approximation less appropriate. Other studies e.g., have also shown that aerosols and clouds further limit the accuracy of the geometrical approximation, yet their effect is not investigated in this study. However, in both cases, the VCDs that are calculated with the geometrical approach can be generally considered a relatively good estimation of the VCDs that are obtained by the vertical profiles analyzed in this study. The calculation of VCDs using model-derived dAMFs would possibly further improve the comparison with the profile-derived VCDs. Unfortunately, the dAMFs that were used in do not include longer wavelengths that are appropriate for the dSCDs derived by the MAX-DOAS system used in this study.

Figure B1

Comparison of NO${}_{\mathrm{2}}$ (a) and HCHO (b) VCDs calculated using the geometrical approximation with the VCDs obtained from the integration of the vertical profiles retrieved by MMF (cyan) and MAPA (magenta).

[Figure omitted. See PDF]

Code availability

The MAPA algorithm is available from the authors upon request (Steffen Beirle, [email protected], and Thomas Wagner, [email protected]). The MMF algorithm is also available from the authors upon request (Martina Michaela Friedrich, [email protected]).

Data availability

The datasets that are used in this study are available from the corresponding author upon request.

Author contributions

DK developed the intercomparison and validation strategy of the two inversion algorithms, analyzed the MAX-DOAS data of Thessaloniki, performed the offline retrievals, conducted the data analysis, and wrote the manuscript. MMF provided the MMF source code, supported and guided DK during the whole time for the proper use of the inversion algorithm, and provided a lot of feedback for the interpretation of the OEM-based results. SB and TW provided the MAPA source code along with useful information about the retrievals and contributed to scientific discussions. KAV provided the lidar extinction profiles, IF and AK the Brewer spectrophotometer-derived AOD, and PT the in situ NO${}_{\mathrm{2}}$ data. MVR, FH, and DB reviewed the paper. AB supervised the whole study and provided general guidance for the manuscript preparation. All authors discussed, commented on, and helped review the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Measurement Techniques. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

This study has been partly supported by the FRM4DOAS project (ESA contract no. 4000118181/16/I-EF). We also thank Caroline Fayt ([email protected]) and Thomas Danckaert for the free use of the QDOAS software and Robert Spurr for providing the VLIDORT radiative transfer code package.

Financial support

This research is co-financed by Greece and the European Union (European Social Fund – ESF) through the Operational Programme “Human Resources Development, Education and Lifelong Learning 2014–2020” in the context of the project “Strengthening Human Resources Research Potential via Doctorate Research 2nd Cycle” (MIS 5000432). The research has been also supported by the project “PANhellenic infrastructure for Atmospheric Composition and climatE change” (MIS 5021516), which is implemented under the Action “Reinforcement of the Research and Innovation Infrastructure”, funded by the Operational Programme “Competitiveness, Entrepreneurship and Innovation” (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).

Review statement

This paper was edited by Piet Stammes and reviewed by two anonymous referees.